Numpy: решение системы линейных уравнений методом Гаусса
Пытаюсь в Python создать алгоритм расчета СЛАУ методом Гаусса. Метод заключается в следующем. Составляем матрицу коэффициентов, включая свободный член. Затем приводим матрицу к треугольному виду. Для этого сначала по первому столбцу (с индексом 0) каждый элемент делим на диагональный (a0,0) (в примере — 3,8), вычисляя индекс m, а после пересчитываем строку 2 по формуле: каждый ее элемент (без элемента свободного члена из последнего столбца) минус произведение элемента над ним (из нулевой строки) и индекса m для второй строки. Отдельно отработаем со столбцом свободного члена (здесь алгоритм неважен). Следом аналогичные действия надо проделать для третьей строки элементов (но учитывая, что на первой итерации элементы второй строки преобразованы вышеописанным алгоритмом, а коэффициент m будет считаться по второму столбцу: соответственно делим все его элементы на диагональный элемент из 2-й строки a1,1) (в примере 1,3). Вопрос: я рассчитал вектор-столбец m: m = ([1,000, 1,684, 0,632]) И теперь надо отработать со второй строкой матрицы. И вот здесь сложность с индексацией. Во-первых, не могу перебрать значения m, тип которых float. Во-вторых, неверно задаю индексацию элементов второй строки (по сути — после нулевой это первая строка)
import numpy as np matrix = np.array([[3.8, 6.7, -1.2, 5.2], [6.4, 1.3, -2.7, 3.8], [2.4, -4.5, 3.5, -0.6]]) def gaussFunc(matrix): # расчет len1 (3) и len2 (4) - здесь не приводится # код расчета m по нулевому столбцу: for row in range(len1): for column in range(len2-3): m = matrix[row][column] / matrix[column][column] elem = row-1 # значения столбцов по нулевой строке кладем в # переменную elem for i in range(len(m)-1): # идем циклом по диапазону трех значений m минус #последнее третье: ошибка по len для float while row < (len1-1): # пока строка первая или вторая (в len2 их всего # 3): while column < (len2-1): # пока колонка первая, вторая или третья # (минус столбец свободного # члена): # пересчитанные коэффициенты второй (первой в numpy) строки: # текущий элемент - m по данной строке*верхний элемент в данном # столбце (из строки 0): a = matrix[row][column] - m[i]*matrix[elem][column] 2 ответа 2
В конце приведена ссылка на jupyter notebook с более-менее полным решателем СЛАУ. Плюс трюк, как считать на pyton почти так же быстро, как на Фортране 🙂
Если не обращать внимание на возможное деление на ноль, то привидение к треугольному виду можно записать очень просто:
def gaussFunc(matrix): # функция меняет матрицу через побочные эффекты # если вам нужно сохранить прежнюю матрицу, скопируйте её np.copy for nrow, row in enumerate(matrix): # nrow равен номеру строки # row содержит саму строку матрицы divider = row[nrow] # диагональный элемент # делим на диагональный элемент. row /= divider # теперь надо вычесть приведённую строку из всех нижележащих строчек for lower_row in matrix[nrow+1:]: factor = lower_row[nrow] # элемент строки в колонке nrow lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow # все строки матрицы изменились, в принципе, можно и не возвращать return matrix Результат для вашего примера:
array([[ 1. , 1.76315789, -0.31578947, 1.36842105], [-0. , 1. , 0.06800211, 0.49657354], [ 0. , 0. , 1. , 0.09309401]]) В чём засада. В ходе вычислений может оказаться ноль на диагонали.
matrix = np.array([[1, 0, 0, 1], [0, 0, 1, 2], [0, 1, 0, 3]]) Насколько я помню, перед тем, как делить на диагональный элемент сначала просматривают все строки, начиная с текущей, вниз. Выбирают строку с максимальным значением в текущей колонке и переставляют с текущей. После чего продолжают.
Функция make_identity берёт матрицу, полученную методом Гаусса, и доводит её до единичной. Для этого строки перебираются снизу вверх и вычитаются из вышележащих строк, чтобы обнулить соответствующие колонки.
def make_identity(matrix): # перебор строк в обратном порядке for nrow in range(len(matrix)-1,0,-1): row = matrix[nrow] for upper_row in matrix[:nrow]: factor = upper_row[nrow] upper_row -= factor*row return matrix Результат для вашего примера: make_identity(gaussFunc(np.copy(matrix)))
array([[ 1. , 0. , 0. , 0.53344344], [-0. , 1. , 0. , 0.49024295], [ 0. , 0. , 1. , 0.09309401]]) Вырезаем последний столбец, получим строку корней: roots = make_identity(gaussFunc(np.copy(matrix)))[:,3]
array([0.53344344, 0.49024295, 0.09309401]) Умножаем найденные корни на исходную матрицу и сравниваем с последним столбцом исходной задачи: np.matmul(matrix[. 3], roots.T) - matrix[:,3]
Результат вычисления array([ 0.00000000e+00, -4.44089210e-16, -2.22044605e-16])
Следовательно, корни найдены правильно. С чем и поздравляю.
Метод Гаусса с выбором главного элемента. Плюс добавлена обработка нуля на диагонали.
def gaussPivotFunc(matrix): for nrow in range(len(matrix)): # nrow равен номеру строки # np.argmax возвращает номер строки с максимальным элементом в уменьшенной матрице # которая начинается со строки nrow. Поэтому нужно прибавить nrow к результату pivot = nrow + np.argmax(abs(matrix[nrow:, nrow])) if pivot != nrow: # swap # matrix[nrow], matrix[pivot] = matrix[pivot], matrix[nrow] - не работает. # нужно переставлять строки именно так, как написано ниже matrix[[nrow, pivot]] = matrix[[pivot, nrow]] row = matrix[nrow] divider = row[nrow] # диагональный элемент if abs(divider) < 1e-10: # почти нуль на диагонали. Продолжать не имеет смысла, результат счёта неустойчив raise ValueError(f"Матрица несовместна. Максимальный элемент в столбце : ") # делим на диагональный элемент. row /= divider # теперь надо вычесть приведённую строку из всех нижележащих строчек for lower_row in matrix[nrow+1:]: factor = lower_row[nrow] # элемент строки в колонке nrow lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow # приводим к диагональному виду make_identity(matrix) return matrix В этой функции главный фокус в том, как переставлять строки. Простая формула matrix[nrow], matrix[pivot] = matrix[pivot], matrix[nrow] не работает. При таком присваивании справа стоят указатели на строку, а слева - адреса, куда нужно скопировать значения. То есть при первом присваивании в строку nrow будет скопирована строка pivot , а в строку pivot - содержимое строки nrow -- но оно уже переписано из строки pivot ! То есть фактически строка pivot скопируется в саму себя. И вместо перестановки двух строк будет копия одной строки.
matrix[[nrow, pivot]] = matrix[[pivot, nrow]] - работает. И c явным копированием тоже работает: matrix[nrow], matrix[pivot] = matrix[pivot], np.copy(matrix[nrow])
Помимо собственно решателя дано сравнение с Си-шным решателем numpy.linalg.solve и трюк как ускорить скорость счёта решателя на пайтоне в 20 раз, что почти так же быстро, как чисто Си-шный решатель.
Строго говоря, решатель в numpy даже не Си-шный, а фортрановский LAPACK gesv
Поиск коэффициентов уравнения
В уравнении вида y = alpha + k*x^n необходимо найти коэффициенты alpha , k , n . Показатель степени n от 0 до 1. При использовании функции curve_fit из scipy получаются совсем некорректные значения. Тестовые данные x = [1021.998 10.21998 5.10999] y = [39.8112 15.312 1.5312] . Должно было получится alpha = 215.4611 n = 0.32 k = 4.3831 . если использовать curve_fit, то на выходе alpha = -77523.501 ; k=77518, 284 ; n = 0 .
from scipy.optimize import curve_fit def calc_func(x, alpha, k, n): return alpha + k * x ** n coeff, pcov = curve_fit(calc_func, x[np.isfinite(y)], y[np.isfinite(y)], maxfev=10000) Выяснила, что для такого типа функции такой метод не подходит, нужно решать численными методами минимизации. Есть ли таких случаев метод решения средствами scipy? Может есть другой способ решения
Я бы избавился от степени в уравнении регрессии. Если не ошибаюсь, это получится log(y-alpha) - log(k) - n*log(x) = 0
Можно попробовать. Тогда получается мы линеализируем функцию и возможно метод curve_fit поможет. Спасибо за идею, попробую так.
Добавьте пожалуйста пример данных, на которых получаются некорректные решения. Кроме того, уточните постановку задачи: может ли x быть отрицательным, n - целое, может ли n быть отрицательным, от 0 до 1? Это я к тому, что x^2 и x^(-2) - совсем разные функции, и для них нужно искать приближенное решение по-разному.
1 ответ 1
У вас очень неудачный набор данных. Две точки находятся почти рядом, а третья точка - очень далеко. В результате вычислительная задача получается крайне неустойчивой и curve_fit не сходится.
Для решения можно использовать тот факт, что n меняется в небольшом диапазоне. Тогда вместо решения задачи с тремя параметрами можно итерировать n и вручную выбрать решение.
import numpy as np import scipy.optimize as spo import matplotlib.pyplot as plt x = np.array([1021.998, 10.21998, 5.10999]) y = np.array([39.8112, 15.312, 1.5312]) def fn(x, a,k,n): return a + k*(x**n) for n in np.linspace(0.01,0.3,20): p_1, p_cov_1 = spo.curve_fit(lambda x,a,k: fn(x,a,k, n), x, y ) err = np.sqrt(np.diag(p_cov_1)) print("n=, (a,k)=±, относительная ошибка ".format(n, p_1, err, err/p_1)) Получается вот какой результат:
n=0.010, (a,k)=[-629.213 624.831]±[162.505 156.627], относительная ошибка [-0.258 0.251] n=0.025, (a,k)=[-234.577 230.849]±[64.832 58.94 ], относительная ошибка [-0.276 0.255] n=0.041, (a,k)=[-137.351 134.248]±[40.795 34.886], относительная ошибка [-0.297 0.26 ] n=0.056, (a,k)=[-93.434 90.929]±[29.96 24.033], относительная ошибка [-0.321 0.264] n=0.071, (a,k)=[-68.471 66.538]±[23.823 17.874], относительная ошибка [-0.348 0.269] n=0.086, (a,k)=[-52.407 51.019]±[19.894 13.921], относительная ошибка [-0.38 0.273] n=0.102, (a,k)=[-41.23 40.362]±[17.178 11.18 ], относительная ошибка [-0.417 0.277] n=0.117, (a,k)=[-33.025 32.653]±[15.202 9.175], относительная ошибка [-0.46 0.281] n=0.132, (a,k)=[-26.76 26.862]±[13.708 7.653], относительная ошибка [-0.512 0.285] n=0.147, (a,k)=[-21.834 22.386]±[12.548 6.463], относительная ошибка [-0.575 0.289] n=0.163, (a,k)=[-17.868 18.85 ]±[11.628 5.511], относительная ошибка [-0.651 0.292] n=0.178, (a,k)=[-14.616 16.005]±[10.885 4.737], относительная ошибка [-0.745 0.296] n=0.193, (a,k)=[-11.906 13.684]±[10.278 4.098], относительная ошибка [-0.863 0.299] n=0.208, (a,k)=[-9.621 11.768]±[9.775 3.564], относительная ошибка [-1.016 0.303] n=0.224, (a,k)=[-7.672 10.17 ]±[9.356 3.113], относительная ошибка [-1.219 0.306] n=0.239, (a,k)=[-5.995 8.826]±[9.004 2.729], относительная ошибка [-1.502 0.309] n=0.254, (a,k)=[-4.541 7.688]±[8.706 2.401], относительная ошибка [-1.917 0.312] n=0.269, (a,k)=[-3.27 6.718]±[8.452 2.118], относительная ошибка [-2.585 0.315] n=0.285, (a,k)=[-2.154 5.887]±[8.236 1.873], относительная ошибка [-3.823 0.318] n=0.300, (a,k)=[-1.169 5.171]±[8.05 1.659], относительная ошибка [-6.889 0.321] Лично мне нравится n=0.147 , alpha = -21.834 , k = 22.386 . Абсолютная ошибка и относительная приемлемые. Конечно, не идеальные, но терпимые.
Видите на графике как сильно изгибается кривая? Это от того, что первые две точки находятся практически на одной вертикали, если сравнивать с правой точкой. Для более качественной аппроксимации желательно, чтобы вторая точка была где-то посередине между первой и третьей.
Так как первые две точки почти на одной вертикали, то n должно быть около нуля, чтобы кривая решения быстро росла.
Но в этом случае очень быстро нарастает абсолютная ошибка для alpha и k . Как ни крути, очень неудачные данные.
На хороших данных curve_fit работает вполне себе успешно.
a = -10; k = 3; n = 0.1 u = np.array([1,500,100]) v = fn(u, a,k,n) p, p_cov = spo.curve_fit(fn, u,v) Результат p = array([-10. , 3. , 0.1]) , ковариационная матрица p_cov не определена, так как всего три точки для трёх параметров. Но главное в другом - все три параметра найдены абсолютно точно.