Метод Гаусса — классический метод решения системы линейных алгебраических уравнений (СЛАУ).Рассмотрим систему линейных уравнений с действительными постоянными коэффициентами:
или в матричной форме
Метод Гаусса решения системы линейных уравнений включает в себя 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++
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
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;
}
Результат выполнения
Назад: Алгоритмизация
2
3
4
5
delete[] a[i];
delete[] a;
delete[] x;
delete[] y;
2
// DBL_EPSILON определено во float.h
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.
почему так может быть?
2
3
4
5
6
7
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
где может быть ошибка, давайте вместе разберем?