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

Алгоритмизация / Решение систем линейных уравнений методом Гаусса

 

Метод Гаусса — классический метод решения системы линейных алгебраических уравнений (СЛАУ).Рассмотрим систему линейных уравнений с действительными постоянными коэффициентами:
Система линейных алгебраических уравнений
или в матричной форме
Система линейных алгебраических уравнений в матричной форме
Метод Гаусса решения системы линейных уравнений включает в себя 2 стадии:

  • последовательное (прямое) исключение;
  • обратная подстановка.

Последовательное исключение

Исключения Гаусса основаны на идее последовательного исключения переменных по одной до тех пор, пока не останется только одно уравнение с одной переменной в левой части. Затем это уравнение решается относительно единственной переменной. Таким образом, систему уравнений приводят к треугольной (ступенчатой) форме. Для этого среди элементов первого столбца матрицы выбирают ненулевой (а чаще максимальный) элемент и перемещают его на крайнее верхнее положение перестановкой строк. Затем нормируют все уравнения, разделив его на коэффициент ai1, где i– номер столбца.
Нормирование линейной системы уравнений
Затем вычитают получившуюся после перестановки первую строку из остальных строк:
Нормирование линейной системы уравнений
Получают новую систему уравнений, в которой заменены соответствующие коэффициенты.
Нормирование линейной системы уравнений
После того, как указанные преобразования были совершены, первую строку и первый столбец мысленно вычёркивают и продолжают указанный процесс для всех последующих уравнений пока не останется уравнение с одной неизвестной:
Ступенчатая форма системы линейных уравнений

Обратная подстановка

Обратная подстановка предполагает подстановку полученного на предыдущем шаге значения переменной xn в предыдущие уравнения:
Обратная подстановка
Эта процедура повторяется для всех оставшихся решений:
Обратная подстановка

Иллюстрирующий пример

Пусть дана система уравнений
Система линейных уравнений (пример)
или в матричной форме
Система линейных уравнений в матричной форме (пример)
Выбираем строку с максимальным коэффициентом ai1 и меняем ее с первой.
Перемещение строк в матрице коэффициентов
Нормируем уравнения относительно коэффициента при x1:
Нормирование коэффициентов (пример)
Нормирование уравнений
Вычитаем 1 уравнение из 2 и 3:
Вычитание уравнений
Выбираем строку с наибольшим коэффициентом при ai2 (уравнение 1 не рассматривается) и перемещаем ее на место 2.
Перемещение строк коэффициентов (пример)
Нормируем 2 и 3 уравнения относительно коэффициента при x2


Вычитаем уравнение 2 из 3

Нормируем уравнение 3 относительно коэффициента при x3

Откуда получаем x3=2. Подставляем полученное значение в уравнения 2 и 1 получаем

Подставляя полученное значение x2=5 в уравнение 1, найдем

Таким образом, решением системы уравнений будет вектор
Вектор решения системы линейных уравнений
Реализация на C++

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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include <iostream>
using namespace std;
// Вывод системы уравнений
void sysout(double **a, double *y, int n)
{
  for (int i = 0; i < n; i++) 
  {
    for (int j = 0; j < n; j++) 
    {
      cout << a[i][j] << "*x" << j;
      if (j < n - 1)
        cout << " + ";
    }
    cout << " = " << y[i] << endl;
  }
  return;
}
double * gauss(double **a, double *y, int n) 
{
  double *x, max;
  int k, index;
  const double eps = 0.00001;  // точность
  x = new double[n];
  k = 0;
  while (k < n) 
  {
    // Поиск строки с максимальным a[i][k]
    max = abs(a[k][k]);
    index = k;
    for (int i = k + 1; i < n; i++) 
    {
      if (abs(a[i][k]) > max)
      {
        max = abs(a[i][k]);
        index = i;
      }
    }
    // Перестановка строк
    if (max < eps) 
    {
      // нет ненулевых диагональных элементов
      cout << "Решение получить невозможно из-за нулевого столбца ";
      cout << index << " матрицы A" << endl;
      return 0;
    }
    for (int j = 0; j < n; j++) 
    {
      double temp = a[k][j];
      a[k][j] = a[index][j];
      a[index][j] = temp;
    }
    double temp = y[k];
    y[k] = y[index];
    y[index] = temp;
    // Нормализация уравнений
    for (int i = k; i < n; i++) 
    {
      double temp = a[i][k];
      if (abs(temp) < eps) continue; // для нулевого коэффициента пропустить
      for (int j = 0; j < n; j++) 
        a[i][j] = a[i][j] / temp;
      y[i] = y[i] / temp;
      if (i == k)  continue; // уравнение не вычитать само из себя
      for (int j = 0; j < n; j++)
        a[i][j] = a[i][j] - a[k][j];
      y[i] = y[i] - y[k];
    }
    k++;
  }
  // обратная подстановка
  for (k = n - 1; k >= 0; k--)
  {
    x[k] = y[k];
    for (int i = 0; i < k; i++)
      y[i] = y[i] - a[i][k] * x[k];
  }
  return x;
}
int main() 
{
  double **a, *y, *x;
  int n;
  system("chcp 1251");
  system("cls");
  cout << "Введите количество уравнений: ";
  cin >> n;
  a = new double*[n];
  y = new double[n];
  for (int i = 0; i < n; i++) 
  {
    a[i] = new double[n];
    for (int j = 0; j < n; j++) 
    {
      cout << "a[" << i << "][" << j << "]= ";
      cin >> a[i][j];
    }
  }
  for (int i = 0; i < n; i++) 
  {
    cout << "y[" << i << "]= ";
    cin >> y[i];
  }
  sysout(a, y, n);
  x = gauss(a, y, n);
  for (int i = 0; i < n; i++) 
    cout << "x[" << i << "]=" << x[i] << endl;
  cin.get(); cin.get();
  return 0;
}

Результат выполнения
Пример


Назад: Алгоритмизация

Комментариев к записи: 68

  • Добрый день. Спасибо за реализацию метода. В строках 60 и 64 можно перебирать j не от 0, а от k. так как от 0 до k будут 0. Их нормировать и вычитать бессмысленно. На больших матрицах чуток ускорится работа :D

  • Скажите пожалуйста, есть ли на Вашем сайте реализация метода Крамера (https://www.mathros.net.ua/rozvjazok-systemy-linijnyh-algebraichnyh-rivnjan-metodom-kramera.html)?

    • Елена Вставская
      На сайте пока нет. Но код реализации имеется. Если нужен, пишите в личку. Возможно, размещу на сайте при случае.

  • Александра
    А есть ли какая-то функция, которая будет выводить уже преобразованную систему с нулевыми коэф-ми ниже главной диагонали?

  • Евгений
    Спасибо за объяснение. Стало понятно как реализовать подобный алгоритм.



    • Вячеслав
      Можно не искать и соответственно не производить замену строк основной матрицы и элементов правого вектора, но тогда это будет решением без выбора ведущего элемента по столбцу. Для СЛАУ малой размерности это чаще всего не столь критично, особенно если все коэффициенты системы примерно равны по модулю, но если приходится решать достаточно большие СЛАУ это приводит к резкой потере точности решения (для очень больших СЛАУ алгоритм без выбора ведущего элемента вообще не стоит рассматривать). Очень примерно суть в том, что бы в процессе решения не делить большие числа на малые (при поиске обнуляющих множителей), из-за чего быстрее накапливается вычислительная погрешность, это уже немного другая тема для изучения.

  • Дмитрий
    Почему ведущим коэффициентом a чаще выбирают максимальное значение, а не просто ненулевой?

    • Елена Вставская
      Треугольная матрица, к которой в итоге приводится система, не должна быть вырожденный, иначе она не имеет решения. Поэтому уравнения переставляют в таком порядке, чтобы ведущие коэффициенты были ненулевым.

  • А зачем вообще в принципе проверять условие отличия от 0 диагональных элементов матрицы ?

    • Елена Вставская
      Мы не сможем составить треугольную матрицу, если будут нулевые диагональные элементы, и такая система не решится

    • Какие есть варианты кода на решение пепеопределенных СЛАУ методом Гаусса?

  • Уважаемая Елена, В строке 22 Вы определяете точность: const double eps = 0.00001; А в строке 39 в условии if (max < eps) определяете, что получить решение невозможно, тк все элементы в столбце равны нулю (ну, по Вашему условию). По моему, такое определение точности не имеет смысла: переменные типа double могут хранить значения до 1е-308 по модулю, для этого они занимают 64 бита. Да и 32 битный float хранит значения до порядка 1е-34. Своим условием в строке 39 Вы всю точность и дабла и флота ограничиваете до 1е-5. По этому определению матрица, все элeменты которой уже порядка 1е-6, наполенена нулями, я не говорю про меньшие значения. А расчеты могут вестись на очень малых числах, разное в жизни бывает.

    • Елена Вставская
      Если все значения в матрице очень малые, то можно их промасштабировать, вместе с век гром свободных членов. При необходимости поменять значение точности. Разные есть варианты.

  • Как говорилось в комментариях, это метод Жордана-Гаусса, но при выводе все числа над главной диагональю не нулевые. Есть ли решение данной проблемы? Или просто за долгое время этот код плавно перешел в метод Гаусса?

    • Елена Вставская
      Это метод Гаусса. Поэтому и коэффициенты над главной диагональю ненулевые. Если бы система решалась методом Жордана-Гаусса, то ее необходимо было бы привести к канонической форме, когда коэффициенты на главной диагонали равны 1, а все остальные коэффициенты равны 0.

      • Могли бы дать какие-то подсказки как привести этот код к методу Жордана-Гаусса?

        • Елена Вставская
          Для метода Жордана-Гаусса не будет обратной подстановки (строки 70-76 примера). Вместо этого нужно получить диагональную матрицу - продолжить исключать все элементы, не лежащие на главной диагонали.

  • Валерий
    В программе ошибка. В строке 77 должно быть return y; Не будьте пожирателями времени - исправьте!

    • Елена Вставская
      Пока не увидела проблемы. Мы определяем вектор x, который и возвращает функция.

  • А как проверить систему на совместность и на случай бесконечного множества решений? Помогите, пожалуйста

  • Дмитрий
    Можете подсказать, в цикле while значит переменная k? Как я понял, это номер столбца, это верно?)

  • А вы в курсе, что вы использовали метод Жордана-Гаусса? И что это две разные вещи?

  • Здраствуйте не подскажите как решить уравнение с вырожденной матрицей (нету решения, либо множество) или какой алгоритм позволяет найти решения.

    • Елена Вставская
      С вырожденный матрицей действительно будет множество решений

  • Все сделано замечательно, только небольшое замечание Елена. Ваша программа не будет выполнятся корректно, если система линейных уравнений, имеет бесконечное множество решений - тут Вам нужно подумать что возвращать при if (max < eps). С уважением, Олег

  • 1
    2
    3
    4
    5
    for (int i = 0; i < n; i++)
    delete[] a[i];
    delete[] a;
    delete[] x;
    delete[] y;
    Дабы не было утечек памяти.


    • Елена Вставская
      Да, есть такая недоработка. Проверка совместности системы в примере отсутствует.

  • Вместо
    1
    const double eps = 0.00001;  // точность
    можно прописать
    1
    2
    const double eps = DBL_EPSILON;  // точность
    // DBL_EPSILON определено во float.h

    • Елена Вставская
      Да, но для одной константы придется вставлять в проект целую библиотеку float.h.

      • Это не библиотека, нет идущих вместе никаких исполняемых модулей, это просто определения констант, точнее макросов препроцессора, т.е. вместо вашего 0.00001 просто подставится число, в итоге даже размер исполняемого файла не изменится.


  • Алексей
    пишет [Error] 'abs' was not declared in this scope, объявляю и появляется куча ошибок, можете сказать как поправить?

    • Елена Вставская
      Попробуйте подключить <cmath> или воспользоваться функцией fabs().

  • Дмитрий
    А есть ли какой-то способ решить систему из 2х уровнений, но с 3мя неизвестными ? К примеру, если нам дана сумма их a+b+c=const(s) и const(1)*a+const(2)*b+const(3)*c = const(s)*const(4) ? Ну или, что-то подобное, Спасибо.

    • Дмитрий, для того чтобы система имела единственное решение необходимым условием является одинаковое количество уравнений и неизвестных. В случае если количество неизвестных больше, чем количество уравнений, система будет иметь бесконечное множество решений.

      • Добрый день, а можно ли перепрограммировать текущий код, так чтобы программа могла также решать системы с неодинаковым количеством переменных и количеством строк, как например сделали на этом сайте: http://math.semestr.ru/gauss/gauss.php ? С Уважением, Игорь.

        • Елена Вставская
          Если количество переменных больше количества уравнений, то часть переменных, превышающих количество уравнений, задается в качестве свободных. Если количество уравнений превышает количество переменных, то такая система вообще может не иметь решения в случае если она не является вырожденной. Пример вырожденной системы уравнений: x1 + 2*x2 = 3 3*x1 + 6*x2 = 9

  • Здравствуйте! Спасибо за код. У меня вылазят ошибки sysout: идентификатор не найден и gauss: идентификатор не найден хотя все сделал, как у вас. В чем может быть хитрость?


      • Ответить без Вашего кода довольно сложно. У меня приведенный пример компилируется и работает.

  • Анастасия
    Можно еще вопрос? в моей системе программа не выдает 0. т.е. у меня система:
    10x1-7x2=7;
    -3x1+2x2+6x3=4;
    5x1-x2+5x3=6.
    ответ, посчитанный вручную :
    x1=0;
    x2=-1;
    x3=1.
    программа выдает:
    x1=1.11022e-016;
    x2=-1;
    x3=1.
    почему так может быть?

  • Добрый день, после x = gauss(a, y, n); в главной функции main, нужно добавить проверку исключения, а именно:
    1
    2
    3
    4
    5
    6
    7
    sysout(a, y, n);
    x = gauss(a, y, n);
    if (x != 0)
    {
    for (int i = 0; i < n; i++) {
    cout << "x[" << i << "]=" << x[i] << endl;
    }
    еще меня вот что смутило, например есть такая система уравнений: (в скобочках я имею в виду нижний регистр)
    x(2)+3x(3)-x(4)=2
    2*x(1)-3*x(2) +2*x(3)=-2
    2*x(1)-4*x(2)+x(4)=3
    -2*x(1)+5*x(2)-3*x(3)+3*x(4)=1
    результат вычисления на бумаге у меня получился такой:
    x1=23,3
    x2=-10,2
    x3=7
    x4=8.8
    в программе такой:
    x[0]=-2.25833
    x[1]=-1.25
    x[2]=1.38333
    x[3]=2.8
    где может быть ошибка, давайте вместе разберем?

  • Добрый день, а если бы, например, было 4 уравнения с 5-ю неизвестными (x1,x2,x3,x4,x5), то это можно запрограммировать? (так как вроде по формулировке - бесконечное множество решений) С Уважением, Игорь.

    • Да, так и есть. Количество уравнений должно соответствовать количеству неизвестных, причем система должна быть невырожденной, чтобы она имела единственное решение.

  • Добрый день, почему после нормирования 2 и 3 уравнения относительно коэффициента при x2 в первой строке вместо [1 0,4 0,2] * [x1] = [9,4] стало [1 0,5 0,25]*[x1] = [10]? С Уважением, Игорь.






        • Для однозначного решения количество уравнений должно соответствовать количеству неизвестных. Поэтому для прямоугольных матриц однозначного решения не существует.


          • Все равно не поняла вопроса. Какую задачу требуется решить?

          • Z2    Решение СЛАУ по модулю, по модулю 2 Z2  Z2 

          • Тут метод Гаусса вряд ли поможет. Это - тема для отдельной статьи.

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *