python — Numpy: решение системы линейных уравнений методом Гаусса
В конце приведена ссылка на 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])
Следовательно, корни найдены правильно. С чем и поздравляю.
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"Матрица несовместна. Максимальный элемент в столбце {nrow}: {divider:.3g}") # делим на диагональный элемент. row /= divider # теперь надо вычесть приведённую строку из всех нижележащих строчек for lower_row in matrix[nrow+1:]: factor = lower_row[nrow] # элемент строки в колонке nrow lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow # приводим к диагональному виду make_identity(matrix) return matrix
В этой функции главный фокус в том, как переставлять строки. Простая формула
не работает. При таком присваивании справа стоят указатели на строку, а слева — адреса, куда нужно скопировать значения. То есть при первом присваивании в строку 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
Jupyter Notebook с кодом решателя СЛАУ
Помимо собственно решателя дано сравнение с Си-шным решателем numpy.linalg.solve
и трюк как ускорить скорость счёта решателя на пайтоне в 20 раз, что почти так же быстро, как чисто Си-шный решатель.
Строго говоря, решатель в numpy
даже не Си-шный, а фортрановский LAPACK gesv
НОУ ИНТУИТ | Лекция | Компьютерное моделирование и решение линейных и нелинейных многомерных систем
< Лекция 8 || Лекция 9: 123 || Лекция 10 >
Аннотация: Лекция рассматривает метод и алгоритм решения систем линейных уравнений методом Гаусса
Ключевые слова: оптимальный план, транспортная, гипотеза, математическая модель, системы линейных уравнений, система линейных уравнений, коэффициенты, свободными членами, матричная форма, метод Гаусса, метод прогонки, коэффициентами системы
При моделировании экономических задач, таких как задачи управления и планирования производства, определения оптимального размещения оборудования, оптимального плана производства, оптимального плана перевозок грузов (транспортная задача), распределения кадров и др., может быть положена гипотеза линейного представления реального мира.
Математические модели таких задач представляются линейными уравнениями. Если задача многомерна, то ее математическая модель представляется системой линейных уравнений.
Линейные математические модели также используются в нелинейных системах при условии, если эта нелинейная система условно линеаризирована.
В общем виде система линейных уравнений имеет вид:
где
aij — коэффициенты при неизвестных системы,
bi — свободные члены,
xj — неизвестные системы,
— номер строки,
— номер столбца,
n — порядок системы.
В матричной форме система линейных уравнений имеет вид:где
Численные методы решения систем линейных уравнений (СЛУ) можно разделить на две группы:
- точные или прямые методы,
- приближенные методы.
Приближенные методы реализуют на ЭВМ нахождение корней с заданной точностью и являются итерационными методами.
Точные методы позволяют получить решение системы за конечное число итераций. К точным методам относятся:
- правило Крамера,
- метод Гаусса,
- метод прогонки.
Решение систем линейных уравнений методом Гаусса
В качестве примера возьмем систему 4 порядка.
( 9.1) |
— ведущий элемент первой строки.
Если , то
( 9.2) |
Обозначим:
( 9.3) |
( 9. 4) |
где
Подставляем (9.4) во 2, 3 и 4 уравнение системы (9.1), получим:
Обозначив коэффициенты при неизвестных полученной системы через , а свободные члены через перепишем полученную систему:
( 9.5) |
где
Таким образом, в результате выполнения первого шага прямого хода исходная система (9.1) n-го порядка преобразована к совокупности уравнения (9.4) и системы линейных уравнений (9.5), порядок которой равен n-1.
На втором шаге прямого хода (к=2) из первого уравнения системы (9. 5) находим x2.
-ведущий элемент первой строки системы (9.5).
Если , то из первого уравнения системы (9.5) имеем:
( 9.6) |
где
Подставив выражение (9.6) во второе и третье уравнения системы (9.5), получим новую систему линейных уравнений, порядок которой равен n-2.
( 9.7) |
Таким образом, в результате выполнения второго шага прямого хода исходная система (9.1) преобразована к совокупности уравнений (9.4), (9.6) и системы линейных уравнений (9.7),порядок которой равен n-2.
Дальше >>
< Лекция 8 || Лекция 9: 123 || Лекция 10 >
Пример 1-го решения уравнения диффузии
Пример 1-го решения уравнения диффузииNext: Анализ устойчивости по фон Нейману Up: Уравнение диффузии Предыдущий: Пример одномерной диффузии Теперь решим уравнение диффузии в 1d, используя конечную разность техника обсуждалась выше. Мы ищем решение уравнения. (191) в регионе , при условии первоначального состояние
(200) |
где . Пространственные граничные условия
(201) |
Конечно, мы можем решить эту проблему аналитически дать
(202) |
Обратите внимание, что приведенное выше уравнение описывает гауссову пульс, постепенно уменьшающийся по высоте и расширяется в ширину таким образом, что его площадь сохраняется. Ширина импульса варьируется примерно как
(203) |
Более того, импульс приближается к -функции как .
На рисунке 71 показано сравнение аналитического и численного решений для расчет, выполненный с использованием , , , , и . Видно, что аналитическое и численное решения прекрасно согласуются.
Разумно ожидать, что по мере роста при фиксированном ( т.е. , пространственное разрешение увеличивается при фиксированном временном разрешении) числовое решение должно становиться все более и более точным. Это действительно так, по крайней мере, пока превышает критическое значение. Помимо этого значения существует катастрофический сбой в численном решении. Эта разбивка показана на рис. 72. Видно, что растворе развиваются быстро нарастающие коротковолновые колебания. Действительно, решение в конечном итоге становится фактически бесконечным. Давайте исследуем это необычное и довольно тревожное явление.
Next: Анализ устойчивости по фон Нейману Up: Уравнение диффузии Предыдущий: Пример одномерной диффузии Ричард Фицпатрик 2006-03-29
MathOnWeb — исключение Гаусса
- Что такое система линейных уравнений?
- Некоторые уроки, которые можно извлечь из графического изображения двух уравнений с двумя неизвестными
- Расширенная матрица
- Элементарные операции со строками
- Исключение Гаусса
- Резервный корпус
- Противоречивый случай
Что такое система линейных уравнений?
Линейное уравнение в n неизвестных x 1 , x 2 , … , x n является уравнением вида:
a 1 x 1 + а 2 x 2 + … + a n x n = b
где a 1 , a 90 133 2 , … , a n и b — константы.
Название линейное происходит от того факта, что такое уравнение с двумя неизвестными или переменными представляет собой прямую линию. Система таких уравнений называется система . Пример системы трех линейных уравнений с тремя неизвестными x , y и z это:
Некоторые уроки, которые можно извлечь из построения графика 2 уравнений с 2 неизвестными
Графический метод не очень полезен в качестве вычислительного инструмента, но полезен для визуализации такие понятия, как уникальность решения или значение противоречивых или избыточных систем. Рассмотрим следующую систему двух линейных уравнений с двумя неизвестными:
В этом методе мы просто рисуем графики уравнений, как мы делали справа. Обратите внимание, что график каждого уравнения представляет собой прямую линию. (Это отличительная черта линейной системы. Здесь нет кривых, есть только прямые линии.)
Любая точка на одной прямой является решением одного уравнения, а любая точка на другой прямой является решением другого уравнения. Точка пересечения линий { x = 3,6, y =0,4} является решением которая удовлетворяет обоим уравнениям одновременно. Заметим, что решение единственное. Это потому что линии прямые и есть только одна точка, где они могут пересекаться.
Система линейных уравнений с единственным решением является «нормальной» ситуацией. Однако это можно иметь систему уравнений без решения. Такая система уравнений называется несовместимый . Часто это результат неточного или неправильного анализа физического состояния. система описывается системой уравнений.
Рассмотрим следующую систему двух линейных уравнений с двумя неизвестными:
Эта система уравнений несовместима, поскольку x + y никак не могут равняться 2 и 4 одновременно. На рисунке справа показано, что граф этой системы состоит двух параллельных прямых, которые никогда не пересекаются. Таким образом, решения нет.
Также возможна система уравнений с бесконечным числом решений. Такая система уравнений называется избыточный . Часто это результат неполного анализ физической системы.
Рассмотрим следующую систему уравнений:
Эта система является избыточной, поскольку второе уравнение эквивалентно первому. График состоит из двух линий, лежащих одна над другой. Они «пересекаются» в бесконечном числе точек, поэтому существует бесконечное количество решений.
Подводя итог, линейная система с двумя неизвестными должна иметь как минимум два уравнения, чтобы получить единственное решение. Иметь 1 уравнение недостаточно, потому что 1 уравнение с 2 неизвестными представлено целой строкой. Достаточно двух уравнений, если они не избыточны и не противоречат друг другу. Наличие 3 (или более) уравнений — это слишком много. Третье уравнение должно быть либо избыточным, либо противоречивым.
Эти идеи можно обобщить на линейные системы уравнений с большим количеством неизвестных:
Линейное уравнение с n переменными представляет «гиперплоскость» в n мерном пространстве.
Линейная система уравнений с n неизвестными должна иметь по крайней мере n уравнений, чтобы получить
уникальное решение. Иметь меньшее количество недостаточно; решение не будет единственным.
Достаточно иметь n уравнений, если они не являются избыточными или противоречивыми.
Имея более n уравнений слишком много; система будет либо избыточной, либо непоследовательной.
Расширенная матрица
Мы представим систему уравнений прямоугольным массивом чисел, называемым дополненная матрица . Вот расширенная матрица для приведенного выше примера:
Немного терминологии:
- Элементы расширенной матрицы называются элементами .
- Строки проходят через всю матрицу.
- Столбцы идут вниз по матрице.
- Диагональ матрицы представляет собой набор элементов, который начинается в верхнем, левом углу и идет по диагонали вниз и вправо. Диагональ вышеуказанной матрицы состоит из чисел 4, 1 и 2.
- Любые элементы в позиции a считаются лежащими на выше диагонали , а любой в позиции b ниже диагонали :
Имейте в виду следующее:
- i -й ряд расширенная матрица представляет i -е уравнение
- j -й столбец (слева от вертикальной черты) представляет (коэффициенты) j -я переменная или неизвестная
- вертикальная линия представляет знаки равенства
Элементарные операции с строками
Напомним, что такое уравнение, как:7( x −4)=14,можно решить для x , применив следующие операции:
- Разделив обе части уравнения на одно и то же значение, а именно на 7, получим x -4=2,
- , затем прибавив одинаковое количество к обеим сторонам, а именно 4, чтобы получить х = 6.
Аналогично, решением системы уравнений является любой набор значений всех
переменных, удовлетворяющих всем уравнениям одновременно. Например, система:
имеет решение:
{ х = 7, у = 5, х = 3}.
Это можно проверить, подставив эти значения во все три уравнения и создание трех тождеств.
Система уравнений может быть решена путем обобщения двух операций, описанных выше, и заметив, что решение системы уравнений не изменится, если:
- разделить обе части уравнения на константу, или
- вычитание кратного одного уравнения из другого уравнения.
Эти же операции можно применить к строкам расширенной матрицы, поскольку каждая строка просто представляет уравнение. Затем они называются Elementary Row Operations .
Элементарные операции со строками (E.R.O.):
- E.R.O.#1: Выберите строку расширенной матрицы и разделите (каждый элемент) строку константой.
Пример:
Обозначение означает разделить первую строку расширенной матрицы на 2, чтобы получить новую расширенную матрицу.
- E.R.O.#2: Выберите любую строку расширенной матрицы и вычтите кратное любой другую строку из него (поэлементно).
Пример:
Обозначение означает взять строку 2 и вычесть из нее 3 раза строку 1, чтобы получить новую расширенную матрицу.
Мы будем применять E.R.O. в определенной последовательности (метод исключения Гаусса, описанный ниже) преобразовать расширенную матрицу в треугольную эшелонированную форму . В этой форме расширенная матрица имеет 1 по диагонали, 0 по диагонали и любые числа по диагонали. Например, расширенная матрица:
в виде треугольного эшелона: Эта новая расширенная матрица представляет собой систему уравнений:Она решается обратной подстановкой. Подставляя z = 3 из третьего уравнения в второе уравнение дает y = 5, и подставляя z = 3 и y = 5 в первое уравнение дает x = 7 . Таким образом, полное решение:
{ x = 7, y = 5, z = 3}.
Исключение по Гауссу
В методе исключения по Гауссу элементарные операции со строками (E.R.O.) применяются в определенном чтобы максимально эффективно преобразовать расширенную матрицу в треугольную эшелонированную форму.
Суть метода: Дана система m уравнений в n переменных или неизвестных, выберите первое уравнение и вычтите подходящие множители его из оставшиеся м -1 уравнен. В каждом случае выберите кратное так, чтобы вычитание отменяет или исключает ту же самую переменную, скажем, x 1 . В результате оставшиеся m -1 уравнения содержат только n -1 неизвестных ( x 1 больше не появляется).
Теперь отложите первое уравнение и повторите вышеуказанный процесс с оставшимися м -1 уравнения в n -1 неизвестных.
Продолжайте повторять процесс. Каждый цикл уменьшает количество переменных и количество уравнений. Процесс останавливается, когда:
- Остается одно уравнение с одной переменной. В этом случае существует единственное решение а обратная замена используется для поиска значений других переменных.
- Остались переменные, но нет уравнений. В этом случае нет единственного решения.
- Остались уравнения, но нет переменных (т. е. самые нижние строки расширенной матрицы содержат только нули слева от вертикальной линии). Это свидетельствует о том, что либо система уравнения противоречивы или избыточны. В случае несоответствия информации, содержащейся в уравнениях противоречиво. В случае избыточности все еще может быть уникальное решение и обратная замена может использоваться для поиска значений других переменных.
Примеры всех этих возможностей приведены ниже.
Алгоритм исключения Гаусса:
Преобразование столбцов расширенной матрицы по одному в треугольную ступенчатую форму. Столбец, который в настоящее время преобразуется, называется сводным столбцом . Продолжайте слева направо, пусть основной столбец будет первым столбцом, затем вторым столбцом, и т. д. и, наконец, последний столбец перед вертикальной чертой. Для каждого сводного столбца выполните следующие два шага, прежде чем перейти к следующему сводному столбцу:
- Найдите диагональный элемент в опорном столбце. Этот элемент называется стержнем . Строка, содержащая сводную строку, называется сводной строкой . Разделите каждый элемент в своде ряд по оси (т. е. используйте E.R.O. #1), чтобы получить новую строку оси с 1 в позиции оси.
- Получите 0 в каждой позиции ниже точки поворота, вычитая подходящее кратное значение точки поворота. строку из каждой из строк под ней (т. е. с помощью E.R.O. #2).
По завершении этой процедуры расширенная матрица будет иметь форму треугольного эшелона и можно решить обратной заменой.
Пример: Используйте исключение Гаусса для решения системы уравнений:
Решение: Выполните эту последовательность E.R.O. на расширенной матрице. Установите опорный столбец в столбец 1. Получите 1 в диагональной позиции (подчеркнуто):
Затем установите 0 под опорной точкой (подчеркнуто):
Теперь пусть опорная колонка = вторая колонка. Сначала получите 1 в диагональной позиции:
Затем получите 0 в позиции ниже опорной:
Теперь пусть опорная колонка = третья колонка. Получите 1 в диагональной позиции:
Эта матрица, которая теперь имеет форму треугольного эшелона, представляет:
Она решается обратной подстановкой. Замена z = 3 из третьего уравнения в второе уравнение дает y = 5, и подставляя z = 3 и y = 5 в первое уравнение дает x = 7 . Таким образом, полное решение:
{ x = 7, y = 5, z = 3}.
Пример применения исключения Гаусса к избыточной системе линейных уравнений решите если возможно:
Решение: Выполните эту последовательность E.R.O. на расширенной матрице. Установить сводную колонку в столбец 1. 1 уже находится в опорной позиции, поэтому продолжайте получать 0 под опорной точкой:
Теперь установите опорную колонку на вторую колонку. Сначала получите 1 в диагональной позиции:
Затем получите 0 в позиции ниже опорной:
Теперь установите опорную колонку на третью колонку. Первое, что нужно сделать, это получить 1 в диагональное положение, но нет возможности сделать это. На самом деле эта матрица уже в виде треугольного эшелона и представляет собой:
Эта система уравнений не может быть решена обратной подстановкой, потому что у нас нет значения для z . Последнее уравнение просто утверждает, что 0=0. Единственного решения нет, потому что z могут принимать на любом значении.
Обычно одна или несколько строк нулей внизу расширенной матрицы, которая была помещена в треугольную эшелонированную форму указывает на избыточную систему уравнений.
Пример применения исключения Гаусса к противоречивой системе линейных уравнений
Используйте исключение Гаусса, чтобы представить эту систему уравнений в виде треугольного эшелона решите, если возможно:
Решение: Выполните эту последовательность E.R.O. на расширенной матрице. Установить точку опоры от столбца к столбцу 1. В опорной позиции уже есть 1, поэтому продолжайте получать 0 под опорной точкой:
Теперь установите опорный столбец на второй столбец. В позиции поворота уже есть 1, поэтому продолжайте, чтобы получить 0 ниже опорной:
Теперь установите опорную колонку на третью колонку. Первое, что нужно сделать, это получить 1 в диагональное положение, но нет возможности сделать это. На самом деле эта матрица уже находится в имеет форму треугольного эшелона и представляет собой:
Эта система уравнений противоречива и не имеет решения.