Введение Многие прикладные задачи приводят к необходимости найти общее решение системы нелинейных уравнений.
Общее аналитическое решение системы нелинейных уравнений не найдено.
Есть только численные методы.
Стоит отметить интересный факт: любую систему уравнений над действительными числами можно представить одним эквивалентным уравнением, если принять все уравнения в виде 
, выровняйте их и сложите.
Для численного решения используются итерационные методы последовательных приближений (простая итерация) и метод Ньютона в различных модификациях.
Итерационные процессы естественным образом обобщаются на случай системы нелинейных уравнений вида: 
(1)
Обозначим через 
вектор неизвестных и определим векторную функцию 
Тогда система (1) запишется в виде уравнения: 
(2)
Теперь вернемся к всеми любимому Python и отметим его первенство среди языков программирования, которые люди хотят изучать [1].
Этот факт является дополнительным стимулом рассмотреть численные методы на Python. Однако среди энтузиастов Python бытует мнение, что специальные библиотечные функции, такие как scipy.optimize.root, spsolve_trianular, newton_krylov , являются лучшим выбором для решения задач численными методами.
С этим трудно не согласиться хотя бы потому, что разнообразие модулей подняло Python на вершину популярности.
Однако бывают случаи, когда даже при поверхностном рассмотрении использование прямых известных методов без использования специальных функций библиотеки SciPy также дает хорошие результаты.
Другими словами, новое – это хорошо забытое старое.
Так, в публикации [2] на основе вычислительных экспериментов доказано, что библиотечная функция newton_krylov, предназначенная для решения больших систем нелинейных уравнений, имеет вдвое производительность алгоритма TSLS+WD.
(двухшаговый метод наименьших квадратов), реализованный с использованием библиотеки NumPy. Цель данной публикации представляет собой сравнение по количеству итераций, производительности, а главное, результата решения модельной задачи в виде системы ста нелинейных алгебраических уравнений с использованием библиотечной функции scipy.optimize.root и метода Ньютона, реализованного используя библиотеку NumPy.
Возможности решателя scipy.optimize.root для численного решения систем алгебраических нелинейных уравнений
В качестве базы сравнения была выбрана библиотечная функция scipy.optimize.root, поскольку она имеет обширную библиотеку методов, подходящих для сравнительного анализа.scipy.optimize.root ( fun, x0, args=(), метод='hybr', jac=None, tol=None, обратный вызов=None, ptions=None ) веселье — Векторная функция для поиска корня.
х0 –Начальные условия для поиска корней метод: гибрид - используется модификация гибридного метода Пауэлла; лм – решает системы нелинейных уравнений методом наименьших квадратов.
Как следует из документации [3] методы Бройден1, Бройден2, Андерсон, линейное смешивание, диагбройден, захватывающее смешивание, Крылов являются точными методами Ньютона.
Остальные параметры являются «необязательными» и их можно найти в документации.
Методы решения систем нелинейных уравнений
Следующий материал действительно можно прочитать в литературе, например в [4], но я уважаю своего читателя и для его удобства изложу вывод метода в сокращенном виде, если это возможно.Те, кто не любит формулы , этот раздел пропускается.
В методе Ньютона новое приближение решения системы уравнений (2) определяется из решения системы линейных уравнений : 
(3)
Определим матрицу Якобиана: 
(4)
Запишем (3) в виде: 
(5)
Многие одношаговые методы приближенного решения (2) по аналогии с двухслойными итерационными методами решения систем линейных алгебраических уравнений можно записать в виде: 
(6)
Где 
— параметры итерации, a 
представляет собой квадратную матрицу размера n x n с обратной.
При использовании обозначений (6) метод Ньютона (5) соответствует выбору: 
Система линейных уравнений (5) для нахождения нового приближения 
можно решить итеративно.
В данном случае мы имеем двухэтапный итерационный процесс с внешними и внутренними итерациями.
Например, внешний итерационный процесс может осуществляться методом Ньютона, а внутренние итерации – на основе итерационного метода Зейделя При решении систем нелинейных уравнений можно использовать прямые аналоги стандартных итерационных методов, которые применяются для решения систем линейных уравнений.
Нелинейный метод Зейделя, примененный к решению (2), дает: 
(7)
При этом каждая компонента нового приближения из решения нелинейного уравнения может быть получена на основе простого итерационного метода и метода Ньютона в различных модификациях.
Таким образом, мы вновь приходим к двухэтапному итерационному методу, в котором внешние итерации осуществляются в соответствии с методом Зейделя, а внутренние итерации – в соответствии с методом Ньютона.
Основные вычислительные трудности использования метода Ньютона для приближенного решения систем нелинейных уравнений связаны с необходимостью решения линейной системы уравнений с матрицей Якоби на каждой итерации , и эта матрица меняется от итерации к итерации.
В модифицированном методе Ньютона матрица Якоби инвертируется только один раз: 
(8)
Выбор функции модели
Подобный выбор является непростой задачей, так как при увеличении числа уравнений в системе в соответствии с увеличением числа переменных результат решения не должен меняться, так как в противном случае невозможно проследить правильность решения задачи.система уравнений при сравнении двух методов.
Я даю следующее решение для модельной функции:
Функция f создает систему n нелинейных уравнений, решение которой не зависит от количества уравнений и равно одному для каждой из n переменных.n=100 def f(x): f = zeros([n]) for i in arange(0,n-1,1): f[i] = (3 + 2*x[i])*x[i] - x[i-1] - 2*x[i+1] - 2 f [0] = (3 + 2*x[0] )*x[0] - 2*x[1] - 3 f[n-1] = (3 + 2*x[n-1] )*x[n-1] - x[n-2] - 4 return f
Программа для тестирования на модельной функции с результатами решения системы алгебраических нелинейных уравнений с использованием библиотечной функцииоптимизировать.
root для разных способов поиска корней from numpy import*
from scipy import optimize
import time
ti = time.clock()
n=100
def f(x):
f = zeros([n])
for i in arange(0,n-1,1):
f[i] = (3 + 2*x[i])*x[i] - x[i-1] - 2*x[i+1] - 2
f [0] = (3 + 2*x[0] )*x[0] - 2*x[1] - 3
f[n-1] = (3 + 2*x[n-1] )*x[n-1] - x[n-2] - 4
return f
x0 =zeros([n])
sol = optimize.root(f,x0, method='krylov')
print('Solution:\n', sol.x)
print('Krylov method iteration = ',sol.nit)
print('Optimize root time', round(time.clock()-ti,3), 'seconds')
Прошел только один из методов, приведенных в документации [3].
тестирование по результату решения модельной функции, это метод «Крылова» .
Решение для n=100: Решение: [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.] Итерация метода Крылова = 4219 Оптимизировать время root 7,239 секунды: Решение для n=200 Решение: [1.00000018 0.99999972 0.99999985 1.00000001 0.99999992 1.00000049 0.99999998 0.99999992 0.99999991 1.00000001 1.00000013 1.00000002 0.9999997 0.99999987 1.00000005 0.99999978 1.0000002 1.00000012 1.00000023 1.00000017 0.99999979 1.00000012 1.00000026 0.99999987 1.00000014 0.99999979 0.99999988 1.00000046 1.00000064 1.00000007 1.00000049 1.00000005 1.00000032 1.00000031 1.00000028 0.99999992 1.0000003 1.0000001 0.99999971 1.00000023 1.00000039 1.0000003 1.00000013 0.9999999 0.99999993 0.99999996 1.00000008 1.00000016 1.00000034 1.00000004 0.99999993 0.99999987 0.99999969 0.99999985 0.99999981 1.00000051 1.0000004 1.00000035 0.9999998 1.00000065 1.00000061 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.0000006 1.00000059 1.00000056 1.00000047 1.00000016 1.00000018 0.99999988 1.00000061 1.00000002 1.00000033 1.00000034 1.0000004 1.00000046 1.00000009 1.00000024 1.00000017 1.00000014 1.00000054 1.00000006 0.99999964 0.99999968 1.00000005 1.00000049 1.0000005 1.00000028 1.00000029 1.00000027 1.00000027 0.9999998 1.00000005 0.99999974 0.99999978 0.99999988 1.00000015 1.00000007 1.00000005 0.99999973 1.00000006 0.99999995 1.00000021 1.00000031 1.00000058 1.00000023 1.00000023 1.00000044 0.99999985 0.99999948 0.99999977 0.99999991 0.99999974 0.99999978 0.99999983 1.0000002 1.00000016 1.00000008 1.00000013 1.00000007 0.99999989 0.99999959 1.00000029 1.0000003 0.99999972 1.00000003 0.99999967 0.99999977 1.00000017 1.00000005 1.00000029 1.00000034 0.99999997 0.99999989 0.99999945 0.99999985 0.99999994 0.99999972 1.00000029 1.00000016] Итерация метода Крылова = 9178 Оптимизировать время root 23,397 секунды Заключение: По мере увеличения количества уравнений вдвое становятся заметными ошибки в решении.
При дальнейшем увеличении n решение становится неприемлемым, что возможно из-за автоматической адаптации к шагу, та же самая причина резкого падения производительности.
Но это только мое предположение.
Программа для тестирования на модельной функции с результатами решения системы алгебраических нелинейных уравнений с помощью программы, написанной на Python 3, с учетом соотношений (1)-(8) для поиска корней модифицированным методом Ньютона
Программа для поиска корней модифицированным методом Ньютона from numpy import*
import time
ti = time.clock()
def jacobian(f, x):
h = 1.0e-4
n = len(x)
Jac = zeros([n,n])
f0 = f(x)
for i in arange(0,n,1):
tt = x[i]
x[i] = tt + h
f1= f(x)
x[i] = tt
Jac [:,i] = (f1 - f0)/h
return Jac, f0
def newton(f, x, tol=1.0e-9):
iterMax = 50
for i in range(iterMax):
Jac, fO = jacobian(f, x)
if sqrt(dot(fO, fO) / len(x)) < tol:
return x, i
dx = linalg.solve(Jac, fO)
x = x - dx
print ("Too many iterations for the Newton method")
n=100
def f(x):
f = zeros([n])
for i in arange(0,n-1,1):
f[i] = (3 + 2*x[i])*x[i] - x[i-1] - 2*x[i+1] - 2
f [0] = (3 + 2*x[0] )*x[0] - 2*x[1] - 3
f[n-1] = (3 + 2*x[n-1] )*x[n-1] - x[n-2] - 4
return f
x0 =zeros([n])
x, iter = newton(f, x0)
print ('Solution:\n', x)
print ('Newton iteration = ', iter)
print('Newton method time', round(time.clock()-ti,3), 'seconds')Решение для n=100: Решение: [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.] Итерация Ньютона = 13 Время метода Ньютона 0,496 секунды Решение для n=200: Решение: [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.] Итерация Ньютона = 14 Время метода Ньютона 1,869 секунды Чтобы убедиться, что программа действительно решает систему, перепишем модельную функцию выхода из корня со значением 1 в виде:
n=10
def f(x):
f = zeros([n])
for i in arange(0,n-1,1):
f[i] = (3 + 2*x[i])*x[i]*sin([i]) - x[i-1] - 2*x[i+1] - 2+e**-x[i]
f [0] = (3 + 2*x[0] )*x[0] - 2*x[1] - 3
f[n-1] = (3 + 2*x[n-1] )*x[n-1] - x[n-2] - 4
return fМы получаем: Решение: [ 0.96472166 0.87777036 0.48175823 -0.26190496 -0.63693762 0.49232062 -1.31649896 0.6865098 0.89609091 0.98509235] Итерация Ньютона = 16 Время метода Ньютона 0,046 секунды Заключение: Программа также работает при изменении функции модели.
Теперь вернемся к исходной функции модели и проверим более широкий диапазон n, например 2 и 500. п=2 Решение: [1. 1.] Итерация Ньютона = 6 Время метода Ньютона 0,048 секунды п=500 п=500 Решение: [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. 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. 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.] Итерация Ньютона = 15 Время метода Ньютона 11,754 секунды.
Выводы:
Программа, написанная на Python с использованием модифицированного метода Ньютона, при решении систем нелинейных уравнений по заданной модельной функции имеет большую устойчивость решения, чем при решении с использованием библиотечной функцииоптимизировать.root(f,x0,method='krylov') для Метод Крылова.
Окончательный вывод относительно скорости работы сделать невозможно из-за разных подходов к ступенчатому управлению.
Ссылки:
- Рейтинг языков программирования 2018.
- Бондарь И.
В.
, Фалейчик Б.
В.
Безматричные итерационные процессы с подавлением среднеквадратической ошибки для больших систем нелинейных уравнений.
- scipy.optimize.root.
- Вабищевич П.
Н.
Численные методы: Вычислительный практикум.
– М.
: Книжный дом «ЛИБРОКОМ», 2010. – 320 с.
-
Как Я Писал Интернет-Магазин В 15 Лет
19 Oct, 24 -
Структура Технических Характеристик
19 Oct, 24 -
Лайфхаки По Копирайтингу
19 Oct, 24
