Пытаюсь в 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]
Отслеживать
12.4k 1 1 золотой знак 6 6 серебряных знаков 29 29 бронзовых знаков
задан 17 апр 2021 в 11:08
Alex_Kazantsev Alex_Kazantsev
569 1 1 золотой знак 5 5 серебряных знаков 18 18 бронзовых знаков
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
Результат для вашего примера:
Решение слау методом итераций. Метод простых итераций c++.
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])
Следовательно, корни найдены правильно. С чем и поздравляю.
UPDATE
Метод Гаусса с выбором главного элемента. Плюс добавлена обработка нуля на диагонали.
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])
UPDATE 2
Помимо собственно решателя дано сравнение с Си-шным решателем numpy.linalg.solve и трюк как ускорить скорость счёта решателя на пайтоне в 20 раз, что почти так же быстро, как чисто Си-шный решатель.
Строго говоря, решатель в numpy даже не Си-шный, а фортрановский LAPACK gesv
Источник: ru.stackoverflow.com
Лекция 7. Решение задач линейной алгебры с использованием динамических матриц и функций
Рассмотрим решение следующих задач линейной алгебры: ∙ решение систем линейных алгебраических уравнений; ∙ вычисление обратной матрицы; ∙ вычисление определителя.
7.1. Решение систем линейных алгебраических уравнений методом Гаусса
ЗАДАЧА 7.1. Решить систему линейных алгебраических уравнений (7.1).
| a 00 x 0 a 01 x 1 . a 0n − 1 x n − 1 = b 0 | |
| a 10 x 0 a 11 x 1 . a 1n − 1 x n − 1 = b 1 | (7.1) |
| . | |
| a n − 10 x 0 a n − 11 x 1 . a n − 1n − 1 x n − 1 = b n − 1 |
Через матрицу А обозначим матрицу коэффициентов системы, через вектор В – вектор свободных членов, Х – вектор неизвестных:
| a 00 | a 01 | . | a 0n | − 1 | b 0 | x 0 |
| a 10 | a 11 | . | a 1n | − 1 | b 1 | x 1 |
| A = | a n − 10 | . | b = | . | x = | . |
| a n − 11 | . | a n − 1n − 1 | b n − 1 | x n − 1 | ||
Ax = b . Первый этап — прямой ход метода Гаусса , заключается в приведении расширенной матрицы (7.2) к треугольному виду (7.3). Это означает, что все элементы матрицы (7.2) ниже главной диагонали должны быть равны нулю.
| a 00 | a 01 | . a 0n − 1 | b 0 | |||
| A’ = | a 10 | a 11 | . | a 1n − 1 | b 1 | (7.2) |
| a n − 10 | . | a n − 1n − 1 b n − 1 | ||||
| a n − 11 | . | |||||
| a 00 | a 01 . | a 0n − 1 | b 0 | |||
| A’ = | a 11 . | a 1n − 1 | b 1 | (7.3) | ||
| . | ||||||
0 0 . a n − 1n − 1 b n − 1 Для формирования нулевого столбца матрицы (7.3) необходимо из каждой строки (начиная со первой) вычесть нулевую, умноженную на некоторое число М .
В общем виде можно записать так: 1-я строка = 1-я строка – М × 0-я строка 2-я строка = 2-я строка – М × 0-я строка … i -я строка = i -я строка – М × 0-я строка … n-1- я строка = n-1 -я строка – М × 0-я строка Преобразование элементов первой строки будет происходить по формулам:
| a 10 = a 10 − Ma 00 | a 11 = a 11 − Ma 01 … | |
| a 1i = a 1i − Ma 0i | a 1n − 1 = a 1n − 1 − Ma 0n − 1 … | |
| b 1 = b 1 − Mb 0 | ||
| a 10 − Ma 00 = 0 | ||
| M = a 10 | ||
| a 00 | ||
Элементы второй строки и коэффициент М можно рассчитать аналогично:
| a 20 = a 20 − Ma 00 | a 21 = a 21 − Ma 01 … |
| a 2i = a 2i − Ma 0i | a 2n − 1 = a 2n − 1 − Ma 0n − 1 … |
| b 2 = b 2 − Mb 0 | |
| a 20 − Ma 00 = 0 |
M = a 20 a 00 Таким образом, преобразование элементов i –й строки будет происходить следующим образом:
| a i0 = a i0 − Ma 00 | a i1 = a i1 − Ma 01 … |
| a i n − 1 = a i n − 1 − Ma 0 n − 1 | b i = b i − Mb 0 |
a i0 − Ma i0 = 0 Коэффициент М для i –й строки будет равен: M = a ik a kk Далее необходимо повторить описанный выше алгоритм для следующих столбцов матрицы (7.2), причем начинать преобразовывать первый столбец со
второго элемента, второй столбец – с третьего и т.д. В результате система уравнений примет вид (7.4). Алгоритм этого процесса изображен на рис. 7.1.
| 1 |
| k=0; k |
| 2 |
| i=k+1; i |
| 3 |
| M=a[i][k]/a[k][k], j=k; j |
| 4 |
| a[i][j]-=M a[k][j] |
| 5 |
| b[i]-=M b[k] |
Рис 7.1. Блок-схема алгоритма преобразования расширенной матрицы к треугольному виду
| a 00 x 0 a 01 x 1 a 02 x 2 . a 0n − 1 x n − 1 = b 0 | |
| a 11 x 1 a 02 x 2 . a 1n − 1 x n − 1 = b 1 | (7.4) |
| a 22 x 2 . a 1n − 1 x n − 1 = b 1 |
. a n − 1n − 1 x n − 1 = b n − 1
Если в матрице (7.2) на главной диагонали встретится элемент a kk , равный нулю, то расчет коэффициента М для к -й строки будет невозможен. Избежать деления на ноль можно, избавившись от нулевых элементов на главной диагонали.
Для этого перед обнулением элементов в k –м столбце необходимо найти в нем максимальный по модулю элемент (среди расположенных ниже a kk ), запомнить номер строки, в которой он находится, и поменять ее местами с k -й строкой. Алгоритм, отображающий эти преобразования, приведен на рис. 7.2. Рис 7.2. Блок-схема алгоритма перестановки строк расширенной матрицы Решение системы (7.4) называют обратным ходом метода Гаусса . Последнее ( n-1) -е уравнение системы (7.4) имеет вид:
| a n − 1n − 1 x n − 1 = b n − 1 . | |||||||||
| Если | a | ≠ 0 | , то | x | n − 1 | = | b n − 1 | . | |
| n − 1n − 1 | a n − 1n − 1 | ||||||||
| В | случае, если | a n − 1n − 1 = 0 и b n − 1 = 0 , система имеет бесконечное | |||||||
| множество решений. | |||||||||
| При | a n − 1n − 1 = 0 | и | b n − 1 ≠ 0 система решений не имеет. | ||||||
Предпоследнее ( n–2 )-е уравнение системы (7.4) имеет вид a n − 2n − 2 x n − 2 a n − 2n − 1 x n − 1 = b n − 2 . Откуда
x = b n − 2 − a n − 2n − 1 x n − 1 . n − 2 a n − 2n − 2 Следующее ( n –3) -е уравнение системы (7.4) будет выглядеть так: a n − 3n − 3 x n − 3 a n − 3n − 2 x n − 2 a n − 3n − 1 x n − 1 = b n − 3 , и откуда находим:
| x n − 3 = | b n − 3 − a n − 3n − 2 x n − 2 − a n − 3n − 1 x n − 1 | ||||||||
| a n − 3n − 3 | |||||||||
| Отсюда имеем: | |||||||||
| n − 1 | |||||||||
| a | x | a | x | b n − 3 − | ∑ a n − 3 j x j | . | |||
| x n − 3 = b n − 3 − | n − 3n − 2 | n − 2 | n − 2n − 1 | n − 1 | = | j = n − 2 | |||
| a n − 3n − 3 | a n − 3n − 3 | ||||||||
Таким образом, формула для вычисления i -го значения x будет иметь вид:
| n − 1 | |
| b i − ∑ a ij x j | i=n-1. 0 |
| x i = | j = i 1 |
| a ii |
Алгоритм, реализующий обратный ход метода Гаусса, представлен в виде блок-схемы на рис. 7.3.
| 1 | |
| i=n-1; i>=0; i— | |
| 2 | |
| s=0, j=i+1; i | |
| 3 | s+=a[i][j]x[j] |
| 4 | x i =(b i -s)/a ii |
Рис 7.3. Блок-схема обратного хода метода Гаусса На рис.7.4 изображена общая блок-схема метода Гаусса.
| 1 |
| Начало |
| 2 n |
| 3 |
| i=0; i |
| 4 |
| j=0; j |
| 5 |
| a[i][j] |
| 6 |
| b[i] |
| + |
| 20 b[n]=0 |
| + | |
| — | 21 |
| бесконечное | |
| множество | |
| решений | |
| 22 | |
| Нет | |
| решений | |
| 28 | Конец |
| 7 | |
| k=0; k | |
| 8 | |
| max=|a[k][k]| | |
| r=k | |
| 9 | |
| i=k+1; i | |
| 10 | |
| — | |a[i][k]|>max |
| + | |
| 11 | |
| max=|a[i][k]| | |
| 19 | r=i |
| a[n-1][n-1]=0 | |
| — | 12 |
| j=0; j | |
| 23 | 13 |
| i=n-1;i>=0;i— | |
| c=a[k][j] | |
| a[k][j]=a[r][j] | |
| 24 | a[r][j]=c |
| s=0, j=i+1; j | |
| 25 | 14 |
| c=b[k] | |
| s+=a[i][j] x[j] | |
| b[k]=b[r] | |
| b[r]=c | |
| 15 | |
| 26 | i=k+1; i |
| x[i]=(b[i]-s)/a[i][i] | 16 |
| M=a[i][k]/a[k][k],j=k,j | |
| 27 | |
| x[i] | 17 |
| a[i][j]-=M a[k][j] | |
| 18 | |
| b[i]-=M b[k] |
Рис 7.4. Блок-схема метода Гаусса Теперь алгоритм решения СЛАУ, представленный на рис. 7.4 разобьем на главную функцию main() и функцию решения СЛАУ методом Гаусса. В функции main() вводятся исходные данные, затем вызывается функция SLAU решения системы линейных алгебраических уравнений и выводится вектор решения. Функция SLAU предназначена для решения системы линейных алгебраических уравнений методом Гаусса.
При написании функции следует учитывать следующее: в методе Гаусса изменяются матрица коэффициентов и вектор правых частей. Поэтому, для того чтобы их не испортить, в функции SLAU матрицу коэффициентов и вектор правых частей необходимо скопировать во внутренние (рабочие) переменные, и в функции обрабатывать внутренние переменные-копии. Функция SLAU возвращает значение 0, если решение найдено, 1 – если система имеет бесконечное множество решений, 2 – если система не имеет решений. int SLAU(double **matrica_a,int n,double *massiv_b, double *x) < int i,j,k,r; double c,M,max,s, **a, *b; a=new double *[n]; for(i=0;imax) < max=fabs(a[i][k]); r=i; >for(j=0;j c=b[k];b[k]=b[r];b[r]=c; for(i=k+1;i > if (a[n-1][n-1]==0)
if(b[n-1]==0) return -1; else return -2; else < for(i=n-1;i>=0;i—) < for(s=0,j=i+1;j
| > | |
| return 0;> | |
| for(i=0;i | |
| delete [] a[i]; | |
| delete [] a; | |
| delete [] b; | |
| > | |
| int main() | |
| double **a, *b, *x; | cin>>N; |
| cout a=new double *[N]; |
for(i=0;i>a[i][j]; cout>b[i]; result=SLAU(a,N,b,x); if (result==0) < coutelse if (result==-1) cout<<«Great number of Solution»; else if (result==-2) cout<<«No solution»; for(i=0;i
Источник: studfile.net
VIII Международная студенческая научная конференция Студенческий научный форум — 2016
Цель работы – разработка алгоритма решения СЛАУ методом Гаусса и программы, использующей данный алгоритм.
Решение системы линейных алгебраических уравнений (СЛАУ) имеет большое значение, поскольку к нему сводится решение широкого круга сложных практических задач. Особая необходимость в решении СЛАУ возникает при использовании широкого класса моделей и подходов, применяемых при автоматизированном проектировании аппаратуры с учетом электромагнитной совместимости. Так как решение СЛАУ «вручную» представляет собой трудоёмкую задачу, а использование средств математического моделирования (MathCAD, MathLab) требует наличие у пользователя некоторых специальных знаний, представляется актуальной задача создания алгоритма решения СЛАУ и программы, построенной на данном алгоритме.
Под СЛАУ подразумевают систему, содержащую m уравнений и n неизвестных (x1,x2,…,xn).
Параметры aij называют коэффициентами, а bi– свободными членами СЛАУ. Иногда, чтобы подчеркнуть количество уравнений и неизвестных, говорят так «m×n система линейных уравнений», – тем самым указывая, что СЛАУ содержит m уравнений и n неизвестных.
Если все свободные члены равны 0, то СЛАУ называют однородной. Если среди свободных членов есть хотя бы один, отличный от нуля, СЛАУ называют неоднородной.
Решением СЛАУ называют всякую упорядоченную совокупность чисел (α1,α2,…,αn), если элементы этой совокупности, подставленные в заданном порядке вместо неизвестных x1,x2,…,xn, обращают каждое уравнение СЛАУ в тождество.
Если СЛАУ имеет хотя бы одно решение, ее называют совместной, если же решений нет – несовместной. Если совместная СЛАУ имеет ровно одно решение, её именуют определённой, если бесконечное множество решений – неопределённой.
В данной работе, была создана программа на языке C# для решения СЛАУ методом Гаусса.
Словесно, алгоритм работы программы можно представить в следующем виде:
- Формируем массив, в который записывается расширенная матрица из коэффициентов и свободных членов;
- Делим каждый элемент первой строки на коэффициент в первом столбце;
- Вычитаем из всех строк, начиная со второй, первую строку, умноженную на коэффициент в первом столбце;
- Переходим к следующей строке и выполняем пункты 2 и 3, увеличивая номер столбца, из которого мы берем первый элемент;
- Выполняем пункт 4 до тех пор, пока не кончатся строки;
- Последний столбец получившейся матрицы представляет собой столбец корней СЛАУ, притом номер корня равен номеру строки, в котором он находится в преобразованной матрице.
Программный код на языке C#, реализующая решение указанный алгоритм:
public partial class Form1 : Form
Источник: scienceforum.ru