Решение систем линейных уравнений методом жордана гаусса python

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

Источник

Метод Жордана-Гаусса в Python

Метод Жордана-Гаусса, также известный как метод Гаусса-Жордана, является численным методом решения системы линейных уравнений. Он основан на преобразовании расширенной матрицы системы к ступенчатому виду с помощью элементарных преобразований строк.

Алгоритм метода Жордана-Гаусса

1. Выберите первое уравнение и разделите его на ведущий элемент (первый ненулевой элемент слева).
2. Для всех остальных уравнений вычтите из них первое уравнение, умноженное на соответствующий множитель, чтобы их ведущий элемент стал равным нулю.
3. Повторите шаги 1 и 2 для всех уравнений, кроме первого, используя каждый раз следующий элемент на главной диагонали в качестве ведущего элемента.
4. Когда все уравнения станут приведенными, решение системы будет содержаться в последнем столбце расширенной матрицы.

Реализация на Python

def gauss_jordan(matrix): n = len(matrix) for i in range(n): pivot = matrix[i][i] for j in range(i, n + 1): matrix[i][j] /= pivot for k in range(n): if k != i: factor = matrix[k][i] for j in range(i, n + 1): matrix[k][j] -= factor * matrix[i][j] return [row[-1] for row in matrix] matrix = [ [2, -1, 5, 1, -3], [3, 2, 2, -6, -32], [1, 3, 3, -1, -47], [5, -2, -3, 3, 49], ] solution = gauss_jordan(matrix) print("Решение системы уравнений:", solution) 

Выше приведена реализация метода Жордана-Гаусса на языке Python. Данный код можно использовать для решения систем линейных уравнений с произвольным количеством переменных.

• Получить ссылку
• Facebook
• Twitter
• Pinterest
• Электронная почта
• Другие приложения

Комментарии

Отправить комментарий

Популярные сообщения

Python вывести количество элементов списка

Python: Вывод количества элементов списка В этой статье мы рассмотрим как выводить количество элементов списка с помощью языка программирования Python. Использование функции len() Для определения количества элементов в списке в Python, используйте встроенную функцию len() . my_list = [1, 2, 3, 4, 5] elements_count = len(my_list) print("Количество элементов в списке:", elements_count) Этот код создает список my_list , а затем использует функцию len() для подсчета элементов в списке. Результат будет выведен на экран. Использование цикла for Если вы хотите подсчитать количество элементов списка без использования функции len() , вы можете использовать цикл for . my_list = [1, 2, 3, 4, 5] elements_count = 0 for _ in my_list: elements_count += 1 print("Количество элементов в списке:", elements_count) В этом примере мы инициализируем переменную elements_count значением 0, а затем для каждого элемента в списке увел

Python как перевести число в другую систему счисления

Преобразуйте числа как профессионал! Узнайте, как Python может перевести любое число в любую систему счисления. Даже если вы никогда раньше не сталкивались с программированием, эта статья поможет вам стать экспертом в считывании двоичных, восьмеричных и шестнадцатеричных чисел. Не пропустите возможность раскрыть секреты произвольной системы счисления в Python! Python: Перевод числа в другую систему счисления В языке программирования Python преобразование числа в другую систему счисления может быть выполнено с использованием встроенных функций и методов. Преобразование чисел в двоичную систему Python предоставляет встроенную функцию bin() для преобразования числа в двоичную систему. # Пример преобразования числа в двоичную систему num = 18 binary_num = bin(num) print(binary_num) # Вывод: 0b10010 Преобразование чисел в восьмеричную систему Функция oct() в Python преобразует число в восьмеричную систему. # Пример преобразования числа в восьмеричную систему num = 18

Как сделать шашки на python

Как сделать шашки на Python Как сделать шашки на Python В этой статье мы рассмотрим, как создать простую игру в шашки на Python с использованием библиотеки Pygame. Подготовка Для начала установите библиотеку Pygame, используя следующую команду: pip install pygame Создание доски import pygame pygame.init() WIDTH, HEIGHT = 800, 800 ROWS, COLS = 8, 8 SQUARE_SIZE = WIDTH // COLS WHITE = (255, 255, 255) BLACK = (0, 0, 0) RED = (255, 0, 0) BLUE = (0, 0, 255) def draw_board(win): win.fill(WHITE) for row in range(ROWS): for col in range(row % 2, COLS, 2): pygame.draw.rect(win, BLACK, (row * SQUARE_SIZE, col * SQUARE_SIZE, SQUARE_SIZE, SQUARE_SIZE)) def main(): win = pygame.display.set_mode((WIDTH, HEIGHT)) pygame.display.set_caption("Checkers") clock = pygame.time.Clock() run = True while run: clock.tick(60) for event in pygame.event.get(): if event.ty

Источник