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

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

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

{\displaystyle \left\{{\begin{array}{lcr}a_{11} \cdot x_{1}+a_{12} \cdot x_{2}+ \ldots + a_{1n} \cdot x_{n}& = & y_{1}\\ a_{21} \cdot x_{1} + a_{22} \cdot x_{2} + \ldots +a_{2n} \cdot x_{n} & = & y_{2} \\ \ldots & & \\a_{m1} \cdot x_{1} + a_{m2} \cdot x_{2} + \ldots +a_{mn} \cdot x_{n} & =& y_{m}\\\end{array}}\right.}

или в матричной форме

{\displaystyle {\left({\begin{array}{cccc}a_{11} a_{12} & \ldots & a_{1n}\\ a_{21} a_{22} & \ldots & a_{2n}\\ \ldots & & \\ a_{m1} a_{m2} & \ldots & a_{mn}\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}} \right) = \left({\begin{array}{c}y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}} \right) }}

{\displaystyle{A \cdot X = Y}}

Метод Гаусса решения системы линейных уравнений включает в себя 2 стадии:

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

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

Исключения Гаусса основаны на идее последовательного исключения переменных по одной до тех пор, пока не останется только одно уравнение с одной переменной в левой части.

Затем это уравнение решается относительно единственной переменной.

Таким образом, систему уравнений приводят к треугольной (ступенчатой) форме. Для этого среди элементов первого столбца матрицы выбирают ненулевой (а чаще максимальный) элемент и перемещают его на крайнее верхнее положение перестановкой строк.

Затем нормируют все уравнения, разделив его на коэффициент {\displaystyle{a_{i1}}}, где {\displaystyle{i}} – номер столбца.

{\displaystyle{ \left\{{\begin{array}{lcr} x_1 + \displaystyle{\frac {a_{12}} {a_{11}}} \cdot x_2 + & \ldots + & \displaystyle{\frac {a_{1n}} {a_{11}}} \cdot x_n = & \displaystyle{\frac {y_1} {a_{11}}} \\ x_1 + \displaystyle{\frac {a_{22}} {a_{21}}} \cdot x_2 + & \ldots + & \displaystyle{\frac {a_{2n}} {a_{21}}} \cdot x_n = & \displaystyle{\frac {y_2} {a_{21}}} \\ \ldots \\ x_1 + \displaystyle{\frac {a_{m2}} {a_{m1}}} \cdot x_2 + & \ldots + & \displaystyle{\frac {a_{mn}} {a_{m1}}} \cdot x_n = & \displaystyle{\frac {y_m} {a_{m1}}} \\ \end{array}}\right. }}

Затем вычитают получившуюся после перестановки первую строку из остальных строк:

{\displaystyle{ \left\{{\begin{array}{lcr} x_1 + & \displaystyle{\frac {a_{12}} {a_{11}}} \cdot x_2 + & \ldots + & \displaystyle{\frac {a_{1n}} {a_{11}}} \cdot x_n = & \displaystyle{\frac {y_1} {a_{11}}} \\ 0 + & \displaystyle{\left( \frac {a_{22}} {a_{21}} \ — \ \frac {a_{12}} {a_{11}} \right)} \cdot x_2 + & \ldots + & \displaystyle{\left( \frac {a_{2n}} {a_{21}} \ — \ \frac {a_{1n}} {a_{11}} \right)} \cdot x_n = & \displaystyle{\left( \frac {y_2} {a_{21}} \ — \ \frac {y_{1}} {a_{11}} \right)} \\ \ldots \\ 0 + & \displaystyle{\left( \frac {a_{m2}} {a_{m1}} \ — \ \frac {a_{12}} {a_{11}} \right)} \cdot x_2 + & \ldots + & \displaystyle{\left( \frac {a_{mn}} {a_{m1}} \ — \ \frac {a_{1n}} {a_{11}} \right)} \cdot x_n = & \displaystyle{\left( \frac {y_m} {a_{m1}} \ — \ \frac {y_{1}} {a_{11}} \right)} \\ \end{array}}\right. }}

Получают новую систему уравнений, в которой заменены соответствующие коэффициенты.

{\displaystyle{ \left\{{\begin{array}{lcr} x_1 & + a_{12} \prime \cdot x_2 + & \ldots + & a_{1n} \prime \cdot x_n = & y_1 \prime \\ 0 & + a_{22} \prime \cdot x_2 + & \ldots + & a_{2n} \prime \cdot x_n = & y_2 \prime \\ \ldots \\ 0 & + a_{m2} \prime \cdot x_2 + & \ldots + & a_{mn} \prime \cdot x_n = & y_m \prime \\ \end{array}}\right. }}

После того, как указанные преобразования были совершены, первую строку и первый столбец мысленно вычёркивают и продолжают указанный процесс для всех последующих уравнений пока не останется уравнение с одной неизвестной:

{\displaystyle{ \left\{{\begin{array}{lcr} x_1 & + & a_{12} \prime \cdot x_2 &+ & a_{13} \prime \cdot x_3 &+ & \ldots &+ & a_{1n} \prime \cdot x_n & = & y_1 \prime \\ 0 & + & x_2 &+ & a_{23} \prime \prime \cdot x_3 &+ & \ldots &+ & a_{2n} \prime \prime \cdot x_n & = & y_2 \prime \prime \\ 0 & + & 0 &+ & x_3 &+ & \ldots &+ & a_{3n} \prime \prime \prime \cdot x_n & = & y_3 \prime \prime \prime \\ \ldots \\ 0 & + & 0 &+ & 0 &+ & \ldots &+ & x_n & = & y_m^n \prime \\ \end{array}}\right. }}

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

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

{\displaystyle{x_{n-1} = y_{n-1}^{(n-1) \prime} — a_{(n-1)n}^{(n-1) \prime} \cdot x_n}}
 
{\displaystyle{x_{n-2} = y_{n-2}^{(n-2) \prime} — a_{(n-2)n}^{(n-2) \prime} \cdot x_n — a_{(n-2)(n-1)}^{(n-2) \prime} \cdot x_{n-1}}}
 
{\displaystyle{\ldots}}
 
{\displaystyle{x_{2} = y_{2}^{ \prime \prime} — a_{2n}^{\prime \prime} \cdot x_n — a_{2(n-1)}^{\prime \prime} \cdot x_{(n-1)} — \ldots — a_{23}^{ \prime \prime} \cdot x_{3}}}
 
{\displaystyle{x_{1} = y_{1}^{ \prime} — a_{1n}^{\prime} \cdot x_n — a_{1(n-1)}^{\prime} \cdot x_{(n-1)} — \ldots — a_{13}^{ \prime} \cdot x_{3} — a_{12}^{ \prime} \cdot x_{2}}}

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

Пусть дана система уравнений

{\displaystyle \left\{{\begin{array}{lcr}2 \cdot x_{1}+4 \cdot x_{2}+ 1 \cdot x_{3}& = & 36\\ 5 \cdot x_{1} + 2 \cdot x_{2} + 1 \cdot x_{3} & = & 47 \\ 2 \cdot x_{1} + 3 \cdot x_{2} +4 \cdot x_{3} & =& 37\\\end{array}}\right.}

Или в матричной форме

{\displaystyle {\left({\begin{array}{ccc}2 & 4 & 1\\ 5 & 2 & 1\\ 2 & 3 & 4\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}36 \\ 47 \\ 37 \end{array}} \right) }}

Выбираем строку с максимальным коэффициентом {\displaystyle{a_{i1}}} и меняем ее с первой.

{\displaystyle {\left({\begin{array}{ccc}\red{5} & 2 & 1\\ 2 & 4 & 1\\ 2 & 3 & 4\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}47 \\ 36 \\ 37 \end{array}} \right) }}

Нормируем уравнения относительно коэффициента при {\displaystyle{x_1}}:

{\displaystyle {\left({\begin{array}{ccc}1 & \displaystyle{\frac {2}{5}} & \displaystyle {\frac {1} {5}}\\ \\ 1 & \displaystyle {\frac {4}{2}} & \displaystyle {\frac {1}{2}}\\ \\ 1 & \displaystyle {\frac {3}{2}} & \displaystyle {\frac {4}{2}}\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}\displaystyle {\frac {47}{5}} \\ \\ \displaystyle {\frac {36}{2}} \\ \\ \displaystyle {\frac {37}{2}} \end{array}} \right) }}

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

{\displaystyle {\left({\begin{array}{ccc}1 & 0,4 & 0,2\\ 0 & 1,6 & 0,3\\ 0 & 1,1 & 1,8\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}9,4 \\ 8,6 \\ 9,1 \end{array}} \right) }}

Выбираем строку с наибольшим коэффициентом при {\displaystyle{a_{i2}}} (уравнение 1 не рассматривается) и перемещаем ее на место уравнения 2.

{\displaystyle {\left({\begin{array}{ccc}1 & 0,4 & 0,2\\ 0 & \red{1,6} & 0,3\\ 0 & 1,1 & 1,8\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}9,4 \\ 8,6 \\ 9,1 \end{array}} \right) }}

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

{\displaystyle {\left({\begin{array}{ccc}1 & 0,4 & 0,2\\ 0 & 1 & 0,1875\\ 0 & 1 & 1,636\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}9,4 \\ 5,375 \\ 8,272 \end{array}} \right) }}

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

{\displaystyle {\left({\begin{array}{ccc}1 & 0,4 & 0,2\\ 0 & 1 & 0,1875\\ 0 & 0 & 1,4489\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}9,4 \\ 5,375 \\ 2,897 \end{array}} \right) }}

Нормируем уравнение 3 относительно коэффициента при {\displaystyle{x_3}}

{\displaystyle {\left({\begin{array}{ccc}1 & 0,4 & 0,2\\ 0 & 1 & 0,1875\\ 0 & 0 & 1\end{array}}\right) \cdot \left({\begin{array}{c}x_{1} \\ x_{2} \\ x_{3} \end{array}} \right) = \left({\begin{array}{c}9,4 \\ 5,375 \\ 2 \end{array}} \right) }}

Откуда получаем {\displaystyle{x_3=2}}. Подставляем полученное значение в уравнения 2 и 1, получаем

{\displaystyle{x_2=5,375-0,175 \cdot x_3 = 5,375-0,175 \cdot 2 = 5}}
 
{\displaystyle{x_1= 9,4-0,2 \cdot x_3 — 0,4 \cdot x_2 = 9,4-0,2 \cdot 2 — 0,4 \cdot 5 = 7}}

Таким образом, решением системы уравнений будет вектор

{\displaystyle{X=(7 \ \ 5 \ \ 2)^T}}

Реализация на 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
110
111
112
113
114
#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 = k; 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;
  for (int i = 0; i < n; i++)
    delete[] a[i];
  delete[] a;
  delete[] y;
  delete[] x;
  cin.get(); cin.get();
  return 0;
}

Результат выполнения
Метод Гаусса: пример

34 комментария к “Решение систем линейных уравнений методом Гаусса”

  1. Евгений

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

    1. Вячеслав

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

  2. Дмитрий

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

    1. Елена Вставская

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

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

    1. Елена Вставская

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

    2. Андрей

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

  4. Михаил

    Уважаемая Елена,
    В строке 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. Елена Вставская

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

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

    1. Елена Вставская

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

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

        1. Елена Вставская

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

  6. Матвей

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

  7. Дмитрий

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

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

    1. Елена Вставская

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

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

    1. Елена Вставская

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

  10. Эдуард

    Вместо

    1
    const double eps = 0.00001;  // точность

    можно прописать

    1
    2
    const double eps = DBL_EPSILON;  // точность
    // DBL_EPSILON определено во float.h
    1. Елена Вставская

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

      1. Эдуард

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

  11. Алексей

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

    1. Елена Вставская

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

  12. Дмитрий

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

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

      1. Игорь

        Добрый день,

        а можно ли перепрограммировать текущий код, так чтобы программа могла также решать системы с неодинаковым количеством переменных и количеством строк, как например сделали на этом сайте:

        http://math.semestr.ru/gauss/gauss.php

        ?

        С Уважением,
        Игорь.

        1. Елена Вставская

          Если количество переменных больше количества уравнений, то часть переменных, превышающих количество уравнений, задается в качестве свободных.
          Если количество уравнений превышает количество переменных, то такая система вообще может не иметь решения в случае если она не является вырожденной.

          Пример вырожденной системы уравнений:

          x1 + 2*x2 = 3

          3*x1 + 6*x2 = 9

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

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

Прокрутить вверх