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

Численные методы решения систем нелинейных уравнений

Введение

Многие прикладные задачи приводят к необходимости нахождения общего решения системы нелинейных уравнений. Общего аналитического решения системы нелинейных уравнений не найдено. Существуют лишь численные методы.

Следует отметить интересный факт о том, что любая система уравнений над действительными числами может быть представлена одним равносильным уравнением, если взять все уравнения в форме , возвести их в квадрат и сложить.

Для численного решения применяются итерационные методы последовательных приближений (простой итерации) и метод Ньютона в различных модификациях. Итерационные процессы естественным образом обобщаются на случай системы нелинейных уравнений вида:

(1)

Обозначим через вектор неизвестных и определим вектор-функцию Тогда система (1) записывается в виде уравнения:

(2)

Теперь вернёмся к всеми любимому Python и отметим его первенство среди языков программирования, которые хотят изучать [1].

Этот факт является дополнительным стимулом рассмотрения числительных методов именно на Python. Однако, среди любителей Python бытует мнение, что специальные библиотечные функции, такие как scipy.optimize.root, spsolve_trianular, newton_krylov, являются самым лучшим выбором для решения задач численными методами.

С этим трудно не согласится хотя бы потому, что в том числе и разнообразие модулей подняло Python на вершину популярности. Однако, существуют случаи, когда даже при поверхностном рассмотрении использование прямых известных методов без применения специальных функций библиотеки SciPy тоже дают неплохие результаты. Иными словами, новое- это хорошо забытое старое.

Так, в публикации [2], на основании проведенных вычислительных экспериментов, доказано, что библиотечная функция newton_krylov, предназначенная для решения больших систем нелинейных уравнений, имеет в два раза меньшее быстродействие, чем алгоритм TSLS+WD
(two-step least squares), реализованный средствами библиотеки NumPy.

Целью настоящей публикации является сравнение по числу итераций, быстродействию, а главное, по результату решения модельной задачи в виде системы из ста нелинейных алгебраических уравнений при помощи библиотечной функции scipy.optimize.root и методом Ньютона, реализованного средствами библиотеки NumPy.

Возможности решателя scipy.optimize.root для численного решения систем алгебраических нелинейных уравнений

Библиотечная функция scipy.optimize.root выбрана в качестве базы сравнения, потому что имеет обширную библиотеку методов, пригодных для сравнительного анализа.

scipy.optimize.root(fun, x0, args=(), method=’hybr’, jac=None, tol=None,callback=None, ptions=None)
fun — Векторная функция для поиска корня.
x0 –Начальные условия поиска корней

method:
hybr -используется модификация Пауэлл гибридный метод;
lm – решает системы нелинейных уравнений методом наименьших квадратов.
Как следует из документации [3] методы broyden1, broyden2, anderson, linearmixing, diagbroyden, excitingmixing, krylov являются точными методами Ньютона. Остальные параметры являются «не обязательными» и с ними можно ознакомится в документации.

Методы решения систем нелинейных уравнений

Приведенный далее материал действительно можно прочитать в литературе, например в [4], но я уважаю своего читателя и для его удобства приведу вывод метода по возможности в сокращенном виде. Те, кто не любит формулы, этот раздел пропускают.

В методе Ньютона новое приближение для решения системы уравнений (2) определяется из решения системы линейных уравнений:

(3)

Определим матрицу Якоби:

(4)

Запишем(3) в виде:

(5)

Многие одношаговые методы для приближенного решения (2) по аналогии с двухслойными итерационными методами для решения систем линейных алгебраических уравнений можно записать в виде:

(6)

где — итерационные параметры, a — квадратная матрица n х n, имеющая обратную.

При использовании записи (6) метод Ньютона (5) соответствует выбору:

Система линейных уравнений (5) для нахождения нового приближения может решаться итерационно. В этом случае мы имеем двухступенчатый итерационный процесс с внешними и внутренними итерациями. Например, внешний итерационный процесс может осуществляться по методу Ньютона, а внутренние итерации — на основе итерационного метода Зейделя

При решении систем нелинейных уравнений можно использовать прямые аналоги стандартных итерационных методов, которые применяются для решения систем линейных уравнений. Нелинейный метод Зейделя применительно к решению (2) дает:

(7)

В этом случае каждую компоненту нового приближения из решения нелинейного уравнения, можно получить на основе метода простой итерации и метода Ньютона в различных модификациях. Тем самым снова приходим к двухступенчатому итерационному методу, в котором внешние итерации проводятся в соответствии с методом Зейделя, а внутренние — с методом Ньютона.

Основные вычислительные сложности применения метода Ньютона для приближенного решения систем нелинейных уравнений связаны с необходимостью решения линейной системы уравнений с матрицей Якоби на каждой итерации, причем от итерации к итерации эта матрица меняется. В модифицированном методе Ньютона матрица Якоби обращается только один раз:

(8)

Выбор модельной функции

Такой выбор не является простой задачей, поскольку при увеличении числа уравнений в системе в соответствии с ростом числа переменных результат решения не должен меняться, поскольку в противном случае невозможно отследить правильность решения системы уравнений при сравнении двух методов. Привожу следующее решение для модельной функции:

Функция f создаёт систему из n нелинейных уравнений, решение которой не зависит от числа уравнений и для каждой из n переменных равно единице.

Программа для тестирования на модельной функции c результатами решения системы алгебраических нелинейных уравнений с помощью библиотечной функции optimize.root для разных методов отыскания корней

Только один из методов, приведенных в документации [3] прошёл тестирование по результату решения модельной функции, это метод ‘krylov’.

Решение для n=100:

Вывод: С увеличением числа уравнений вдвое заметно появление ошибок в решении. При дальнейшем увеличении n решение становится не приемлемым, что возможно из-за автоматической адаптации к шагу, эта же причина резкого падения быстродействия. Но это только моё предположение.

Программа для тестирования на модельной функции c результатами решения системы алгебраических нелинейных уравнений с помощью программы написанной на Python 3 с учётом соотношений (1)-(8) для отыскания корней по модифицированному методу Ньютона

Решение для n=100:

Решение для n=200:

Solution:
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1.]
Newton iteration = 14
Newton method time 1.869 seconds

Чтобы убедиться в том, что программа действительно решает систему, перепишем модельную функцию для ухода от корня со значением 1 в виде:

Получим:
Solution:
[ 0.96472166 0.87777036 0.48175823 -0.26190496 -0.63693762 0.49232062
-1.31649896 0.6865098 0.89609091 0.98509235]
Newton iteration = 16
Newton method time 0.046 seconds

Вывод: Программа работает и при изменении модельной функции.

Теперь вернёмся к начальной модельной функции и проверим более широкий диапазон для n, например в 2 и 500.
n=2
Solution:
[1. 1.]
Newton iteration = 6
Newton method time 0.048 seconds
n=500

Численные методы решения нелинейных уравнений. Метод Ньютона для решения уравнений с одной переменной

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

Метод был описан Исааком Ньютоном в рукописи De analysi per aequationes numero terminorum infinitas ( лат .О б анализе уравнениями бесконечных рядов), адресованной в 1669 году Барроу , и в работе De metodis fluxionum et serierum infinitarum ( лат.Метод флюксий и бесконечные ряды) или Geometria analytica ( лат.Аналитическая геометрия) в собраниях трудов Ньютона, которая была написана в 1671 году. Однако описание метода существенно отличалось от его нынешнего изложения: Ньютон применял свой метод исключительно к полиномам. Он вычислял не последовательные приближения x_n , а последовательность полиномов и в результате получал приближённое решение x.

Впервые метод был опубликован в трактате Алгебра Джона Валлиса в 1685 году, по просьбе которого он был кратко описан самим Ньютоном. В 1690 году Джозеф Рафсон опубликовал упрощённое описание в работе Analysis aequationum universalis (лат. Общий анализ уравнений). Рафсон рассматривал метод Ньютона как чисто алгебраический и ограничил его применение полиномами, однако при этом он описал метод на основе последовательных приближений x_n вместо более трудной для понимания последовательности полиномов, использованной Ньютоном.

Наконец, в 1740 году метод Ньютона был описан Томасом Симпсоном как итеративный метод первого порядка решения нелинейных уравнений с использованием производной в том виде, в котором он излагается здесь. В той же публикации Симпсон обобщил метод на случай системы из двух уравнений и отметил, что метод Ньютона также может быть применён для решения задач оптимизации путём нахождения нуля производной или градиента.

В соответствии с данным методом задача поиска корня функции сводится к задаче поиска точки пересечения с осью абсцисс касательной, построенной к графику функции .

Рис.1 . График изменение функции

Проведенная в любой точке касательная линия к графику функции определяется производной данной функции в рассматриваемой точке, которая в свою очередь определяется тангенсом угла α ( ). Точка пересечения касательной с осью абсцисс определяется исходя из следующего соотношения в прямоугольном треугольнике: тангенс угла в прямоугольном треугольнике определяется отношением противолежащего катета к прилежащему катету треугольнику. Таким образом, на каждом шаге строится касательная к графику функции в точке очередного приближения . Точка пересечения касательной с осью Ox будет являться следующей точкой приближения . В соответствии с рассматриваемым методом расчет приближенного значения корня на i -итерации производится по формуле:

Наклон прямой подстраивается на каждом шаге наилучшим образом, однако следует обратить внимание на то, что алгоритм не учитывает кривизну графика и следовательно в процессе расчета остается неизвестно в какую сторону может отклониться график.

Условием окончания итерационного процесса является выполнение следующего условия:

где ˗ допустимая погрешность определения корня.

Метод обладает квадратичной сходимостью. Квадратичная скорость сходимость означает, что число верных знаков в приближённом значении удваивается с каждой итерацией.

Математическое обоснование

Пусть дана вещественная функция , которая определена и непрерывна на рассматриваемом участке. Необходимо найти вещественный корень рассматриваемой функции.

Вывод уравнения основано на методе простых итераций, в соответствии с которым уравнение приводят к эквивалентному уравнению при любой функции . Введем понятие сжимающего отображения, которое определяется соотношением .

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

Производная сжимающего отображения определяется в следующем виде:

Выразим из данного выражение переменную при условии принятого ранее утверждения о том, что при необходимо обеспечить условие . В результате получим выражение для определения переменной :

С учетом этого сжимающая функция прием следующий вид:

Таким образом, алгоритм нахождения численного решения уравнения сводится к итерационной процедуре вычисления:

Алгоритм нахождения корня нелинейного уравнения по методу Ньютона для уравнения с одной переменной

1. Задать начальную точку приближенного значения корня функции , а также погрешность расчета (малое положительное число ) и начальный шаг итерации ( ).

2. Выполнить расчет приближенного значения корня функции в соответствии с формулой:

3. Проверяем приближенное значение корня на предмет заданной точности, в случае:

— если разность двух последовательных приближений станет меньше заданной точности , то итерационный процесс заканчивается.

— если разность двух последовательных приближений не достигает необходимой точности , то необходимо продолжить итерационный процесс и перейти к п.2 рассматриваемого алгоритма.

Пример решения уравнений

по методу Ньютона для уравнения с одной переменной

В качестве примера, рассмотрим решение нелинейного уравнения методом Ньютона для уравнения с одной переменной . Корень необходимо найти с точностью в качестве первого приближения .

Вариант решения нелинейного уравнения в программном комплексе MathCAD представлен на рисунке 3.

Результаты расчетов, а именно динамика изменения приближенного значения корня, а также погрешности расчета от шага итерации представлены в графической форме (см. рис.2).

Рис.2 . Результаты расчета по методу Ньютона для уравнения с одной переменной

Для обеспечения заданной точности при поиске приближенного значения корня уравнения в диапазоне необходимо выполнить 4 итерации. На последнем шаге итерации приближенное значение корня нелинейного уравнения будет определяться значением: .

Рис.3 . Листинг программы в MathCad

Модификации метода Ньютона для уравнения с одной переменной

Существует несколько модификаций метода Ньютона, которые направлены на упрощение вычислительного процесса.

Упрощенный метод Ньютона

В соответствии с методом Ньютона требуется вычислять производную функции f(x) на каждом шаге итерации, что ведет к увеличению вычислительных затрат. Для уменьшения затрат, связанных с вычислением производной на каждом шаге расчета, можно произвести замену производной f’( x_n ) в точке x_n в формуле на производную f’(x₀) в точке x₀. В соответствии с данным методом расчета приближенное значение корня определяется по следующей формуле:

Таким образом, на каждом шаге расчета строятся прямые , которые параллельны касательной к кривой y=f(x) в точке B₀ (см. рис.4). Преимуществом данного метода является то, что производная функции вычисляется один раз.

Разностный метод Ньютона

В соответствии с методом Ньютона требуется вычислять производную функции f(x) на каждом шаге итерации, что не всегда удобно, а иногда практически невозможно. Данный способ позволяет производную функции заменить разностным отношением (приближенным значением):

В результате приближенное значение корня функции f(x) будет определяться выражением разностного метода Ньютона:

Двух шаговый метод Ньютона

В результате приближенное значение корня функции f(x) будет определяться следующим выражением:

Метод секущих является двух шаговым, то есть новое приближение определяется двумя предыдущими итерациями и . В методе необходимо задавать два начальных приближения и . Скорость сходимости метода будет линейной.

Для того, чтобы добавить Ваш комментарий к статье, пожалуйста, зарегистрируйтесь на сайте.

5.1. Приближённое решение систем нелинейных уравнений. Метод Ньютона

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

(5.1)

С действительными левыми частями. Систему (5.1) можно представить в матричном виде

(5.2)

Здесь приняты следующие обозначения:

— вектор аргументов, а — вектор – функция.

Для решения системы (5.2) воспользуемся методом последовательных приближений. Предположим, что найдено Р-ое приближение Xp = (X1(P), X2(P) , . Xn(P)) одного из изолированных корней X = (X1, X2, X3, . Xn) векторного уравнения (5.2). Тогда точный корень уравнения (5.2) можно представить в виде

(5.3)

Где — поправка (погрешность) корня на N – ом шаге.

Подставив выражение (5.3) в (5.2), получим

(5.4)

Предположим, что функция F(X) — непрерывно дифференцируема в некоторой выпуклой области, содержащей X и X(P). Тогда левую часть уравнения (5.4) разложим в ряд Тейлора по степеням малого вектора ε(P), ограничиваясь линейными членами:

, (5.5)

Или в развернутом виде:

(5.6)

Из анализа формул (5.5) и (5.6) следует, что под производной F¢(X) следует понимать матрицу Якоби системы функций F1 , F2, . Fn, относительно переменных X1, X2, X3, . Xn, то есть:

. (5.7)

Выражение (5.7) в краткой записи можно представить:

(5.8)

Выражение (5.6) представляет собой линейную систему относительно поправок (I = 1, 2, . N) с матрицей W(X), поэтому формула (5.5) может быть записана в следующем виде:

(5.9)

Отсюда, предполагая, что матрица W(X(P)) — неособенная, получим:

(5.10)

Теперь, подставив выражение (5.10) в формулу (5.3), окончательно получим:

(5.11)

Таким образом, получили вычислительную формулу (метод Ньютона), где в качестве нулевого приближения X(0) можно взять приближенное (грубое) значение искомого корня.

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

(5.12)

Прежде чем разбирать конкретные шаги по решению системы (5.12), распишем в общем виде якобиан для системы из двух уравнений

Здесь A, B, C, D – функционалы от переменных X1, x2. Нас фактически интересует W-1. Пусть матрица W— неособенная, тогда обратная матрица вычисляется

Теперь вернемся к системе (5.12). Графическое решение этой системы дает две точки пересечения: М1 (1.4; -1.5) и М2 (3.4; 2.2). Зададим начальное приближение:

источники:

http://simenergy.ru/math-analysis/solution-methods/45-method-newton-s

http://matica.org.ua/metodichki-i-knigi-po-matematike/chislennye-metody-iu-ia-katcman/5-1-priblizhennoe-reshenie-sistem-nelineinykh-uravnenii-metod-niutona