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
Метод квадратных корней, не до конца вычисляет решение
Реализовывал метод квадратных корней, следовательно у меня при расчетах получались числа с длинными десятичными дробями. Однако помогите плиз разобраться в чем дело:
Почему у меня иксы получаются слишком сокращенными:
x1 = -6.100
x2 = -2.200
x3 = -6.800
x4 = -0.900
x5 = 0.200
То есть только первая цифра десятичной дроби, тогда как, когда какой-нибудь элемент t_arr вывожу, то выводит полное число, с длинными десятичными дробями. Вот код:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
import numpy as np import cmath as cm #Для комплексных чисел import math as m wide_arr = np.array([[1, 3, -2, 0, -2, 0.5], [3, 4, -5, 1, -3, 5.4], [-2, -5, 3, -2, 2, 5.0], [0, 1, -2, 5, 3, 7.5], [-2, -3, 2, 3, 4, 3.3]]) #Формирование матрицы t элементов t_arr = np.zeros((wide_arr.shape), dtype = complex) t_arr[0][0] = cm.sqrt(wide_arr[0][0]) #t[1][1] #Находим t[1][j] for i in range(wide_arr.shape[1]-1): t_arr[0][i] = wide_arr[0][i]/t_arr[0][0] #t[2][2] t_arr[1][1] = cm.sqrt(wide_arr[1][1] - (t_arr[0][1]**2)) #Находим t[2][j] for j in range(2,5): t_arr[1][j] = (wide_arr[1][j] - (t_arr[0][1])*t_arr[0][j])/t_arr[1][1] #t[3][3] t_arr[2][2] = cm.sqrt(wide_arr[2][2] - (t_arr[0][2]**2 + t_arr[1][2]**2) ) #Находим t[3][j] for j in range(3,5): t_arr[2][j] = (wide_arr[2][j] - (t_arr[0][2]*t_arr[0][j] + t_arr[1][2]*t_arr[1][j]))/t_arr[2][2] #t[4][4] t_arr[3][3] = cm.sqrt(wide_arr[3][3] - (t_arr[0][3]**2 + t_arr[1][3]**2 + t_arr[2][3]**2)) #Находим t[4][5] t_arr[3][4] = (wide_arr[3][4] - (t_arr[0][3]*t_arr[0][4] + t_arr[1][3]*t_arr[1][4] + t_arr[2][3]*t_arr[2][4]))/t_arr[3][3] #t[5][5] t_arr[4][4] = cm.sqrt(wide_arr[4][4] - (t_arr[0][4]**2 + t_arr[1][4]**2 + t_arr[2][4]**2 + t_arr[3][4]**2)) #Находим y[i] #y[1] t_arr[0][5] = wide_arr[0][5]/t_arr[0][0] #y[2] t_arr[1][5] = (wide_arr[1][5] - t_arr[0][1]*t_arr[0][5])/t_arr[1][1] #y[3] t_arr[2][5] = (wide_arr[2][5] - (t_arr[0][2]*t_arr[0][5] + t_arr[1][2]*t_arr[1][5]))/t_arr[2][2] #y[4] t_arr[3][5] = (wide_arr[3][5] - (t_arr[0][3]*t_arr[0][5] + t_arr[1][3]*t_arr[1][5] + t_arr[2][3]*t_arr[2][5]))/t_arr[3][3] #y[5] t_arr[4][5] = (wide_arr[4][5] - (t_arr[0][4]*t_arr[0][5] + t_arr[1][4]*t_arr[1][5] + t_arr[2][4]*t_arr[2][5] + t_arr[3][4]*t_arr[3][5]))/t_arr[4][4] #Находим решение x_arr = np.zeros((wide_arr.shape[0],1), dtype = complex) #x[5] x_arr[4] = t_arr[4][5]/t_arr[4][4] #x[4] x_arr[3] = (t_arr[3][5] - t_arr[3][4]*x_arr[4])/t_arr[3][3] #x[3] x_arr[2] = (t_arr[2][5] - (t_arr[2][3]*x_arr[3] + t_arr[2][4]*x_arr[4]))/t_arr[2][2] #x[2] x_arr[1] = (t_arr[1][5] - (t_arr[1][2]*x_arr[2] + t_arr[1][3]*x_arr[3] + t_arr[1][4]*x_arr[4]))/t_arr[1][1] #x[1] x_arr[0] = (t_arr[0][5] - (t_arr[0][1]*x_arr[1] + t_arr[0][2]*x_arr[2] + t_arr[0][3]*x_arr[3] + t_arr[0][4]*x_arr[4]))/t_arr[0][0] def printMatrix ( matrix ): for i in range ( len(matrix) ): for j in range ( len(matrix[i]) ): print ( "".format(matrix[i][j]), end = " ") print('') print('') print('Матрица t,y элементов') printMatrix(t_arr) print('Решение') for i in range(1,6): print('x<> = '.format( i, x_arr[i-1][0].real))
В коде особо копаться не надо, мне бы точность решения иксов повысить. Столько времени потратил, и тут такое
Добавлено через 19 минут
Поэкспериментировал с кодом, поставил вывод до 18 нулей:
print('Решение') for i in range(1,6): print('x<> = '.format( i, x_arr[i-1][0].real))
Решение x1 = -6.100000000000003197 x2 = -2.200000000000000178 x3 = -6.800000000000002487 x4 = -0.900000000000000910 x5 = 0.200000000000000538
Метод квадратных корней
Словарь квадратных корней
Уважаемые гуру. есть вопрос: # Пользователь вводит в цикле целые положительные числа, # пока не.
Подсчитать сумму квадратных корней
Суть задания в том что нужно подсчитать корни но суть в том что изначально идёт корень 1 потом под.
Вывести на экран таблицу квадратных корней
Вывести на экран таблицу квадратных корней первых десяти целых положительных нечетных числе
Создать функцию с неопределенным количеством аргументов находящую сумму квадратных корней чисел
Создать функцию с переменным количеством находящую сумму квадратных корней чисел. Помогите.
Сообщение было отмечено VictorVAlduin как решение
Решение
Уже лучше, но 10 и 11 строки неправильные. Матрицы нужно умножать на обратную, а не делить, т.к. деление поэлементное.
y=np.linalg.inv(L)*B x=np.linalg.inv(L_TRANSPOSE)*y
[[ 18.73077027] [-149.60533654] [ 310.06279236] [-185.08812279]]
Понял
Ну я следовал строго алгоритму из вышеприведенных формул: но тут правда была только транспонированная и обычная матрица
перепутал местами транспонированную и обычную
Список из квадратных корней чисел
#3 Напишите функцию которая принимает на вход список. # Функция создает из этого списка новый.
Написать три алгоритма решения СЛАУ: Метод прогонки, метод квадратных корней, метод вращений
Начал писать курсовую. Нужно написать три алгоритма решения СЛАУ: прогонки, квадратных корней.
Исследовать метод квадратных корней и метод Холецкого для решения СЛАУ
> Исследовать метод квадратных корней и метод Холецкого для решения СЛАУ при суммировании и.
Метод квадратных корней
Задание. Изготовить и испытать подпрограмму для применения метода квадратных корней для решения.
Метод квадратных корней
Добрый день, получил задание написать программу, которая решает СЛАУ данным методом. Программу я.
Метод квадратных корней
Доброго времени суток. Помогите найти ошибку. Действовал по алгоритму, описанному на картинке.
Метод квадратных корней
Ребята, помогите пожалуйста решить СЛАУ методом квадратных корней с точностью 0.001 на Turbo.