Алгоритм метода Гаусса с выбором главного элемента по столбцам.
1. Для m = 1, 2, …, n – 1 выполним преобразования:
Найдем максимальный по абсолютной величине элемент в m-ом столбце. Пусть это будет элемент aim. Если i ≠ m, то меняем местами i-ую и m-ую строки:
r = aij, aij = amj, amj = r, j = 1, …, n; r = bi, bi = bm, bm = r.
Элементы матрицы и вектора после преобразования на m-ом шаге обозначим 
, причем 
.
Делим m-ое уравнение на ведущий элемент 
:
а затем исключаем переменную xm с помощью m-ого уравнения из i-го,
где i = m + 1, …, n:
2. Вычисляем неизвестные xi в обратном порядке:
; 
Приведенный вариант метода Гаусса дает решение и тогда, когда обычный метод Гаусса перестает работать из-за деления на ноль.
При реализации метода Гаусса на каком-либо языке программирования удобно использовать исходные матрицу a и вектор b для хранения промежуточных результатов преобразований. Приведенные формулы преобразований записываются как операторы присваивания, т.е. имена переменных aij и bj записываются без верхних индексов. Ниже приведены различные примеры программ метода Гаусса.
В другом варианте метода, методе Гаусса с выбором главного элемента по строкам на каждом шаге выбирают наибольший по модулю элемент строки и переставляют столбцы так, чтобы ведущий элемент находился на диагонали. В этом варианте метода Гаусса в уравнениях неизвестные переменные меняются местами и в алгоритме надо заботиться о том, чтобы запомнить эти перестановки.
В общем случае метода Гаусса с выбором главного элемента на шаге m ищется максимальный по модулю элемент во всей матрице коэффициентов, и производится перестановка строк и столбцов так, чтобы этот элемент стал диагональным.
Мы не будем рассматривать реализации указанных вариантов метода Гаусса. Отметим, что последний вариант метода Гаусса с выбором главного элемента считается лучшим.
Общее число операций для решения системы уравнений методом Гаусса составляет приблизительно [1] 2n3/3 + 2n2, при этом большая часть, т.е. 2n3/3 операций приходится на прямой ход.
Пример 3.3. Решить методом Гаусса систему уравнений
4x1 + 2x2 – 4x3 + 2x4 = 12,
2x1 + 9x2 – 4x3 – 4x4 = 14,
3x1 + 2x2 – 8x3 + 2x4 = 20,
4x1 + 2x2 + 5x3 + 7x4 = 10.
Решение.Продемонстрируем решение в программе электронных таблиц.
0. Запишем коэффициенты и правые части системы в диапазоне A1:E4.
1. В диапазоне A5:E8 сформируем расширенную матрицу системы, которая получится после исключения переменной x1 из второго, третьего и четвертого уравнений.
В ячейку A5 вводим формулу =A1/$A1; выделим A5 и протянем маркер заполнения (мышью за правый нижний угол) до E5.
Затем в ячейку A6 вводим формулу =A2-A$5*$A2; выделим A6 и протянем маркер заполнения до E6.
Выделим диапазон A6:E6 и протянем маркером заполнения до E8.
В ячейках A5:E8 появятся следующие значения:
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 |
| 0,0000 | 8,0000 | -2,0000 | -5,0000 | 8,0000 |
| 0,0000 | 0,5000 | -5,0000 | 0,5000 | 11,0000 |
| 0,0000 | 0,0000 | 9,0000 | 5,0000 | -2,0000 |
2. В диапазоне A9:E12 сформируем матрицу системы, которая получится после исключения x2 из третьего и четвертого уравнений.
Присвоим строке A9:E9 значения строки A5:E5, т.е. выделим диапазон A9:E9 и введем формулу =A5:E5 и удерживая нажатыми клавиши Ctrl и Shift, нажимаем Enter.
В ячейку B10 вводим формулу =B6/$B6; выделим B10 и протянем маркер заполнения до E10.
Затем в ячейку B11 вводим формулу =B7-B$10*$B7; выделим B11 и протянем мышью маркер заполнения до E11.
Выделим диапазон B11:E11 и протянем маркер заполнения до E12.
В ячейках A9:E12 появятся следующие значения:
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 |
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 |
| 0,0000 | -4,8750 | 0,8125 | 10,5000 | |
| 0,0000 | 9,0000 | 5,0000 | -2,0000 |
3. В диапазоне A13:E16 сформируем матрицу системы, которая получится после исключения x3 из четвертого уравнения.
Перепишем значения строк A9:E10 в строки A13:E14, для этого выделим диапазон A13:E14 и введем формулу =A9:E10 и удерживая клавиши Ctrl и Shift нажатыми, нажимаем Enter.
В ячейку C15 вводим формулу =C11/$C11; выделим C15 и протянем мышью за правый нижний угол до E15.
В ячейку C16 вводим формулу =C12-C$15*$C12, выделим C16 и протянем мышью за правый нижний угол до E16. В ячейках A13:E16 получим значения:
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 |
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 |
| 1,0000 | -0,1667 | -2,1538 | ||
| 0,0000 | 6,5000 | 17,3846 |
4. Сформируем в диапазоне A17:E20 окончательную матрицу системы, которая имеет треугольный вид с единицами на диагонали.
Выделим диапазон A17:E19 и введем формулу =A13:E15 и, удерживая нажатыми клавиши Ctrl и Shift, нажимаем Enter. В ячейку D20 вводим формулу =D16/$D16, выделим C15 и протянем мышью за правый нижний угол до E20. В диапазоне A17:E20 получим:
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 |
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 |
| 0,0000 | 0,0000 | 1,0000 | -0,1667 | -2,1538 |
| 1,0000 | 2,6746 |
В таблице 3.1 показаны введенные формулы, которые можно увидеть, если выполнить команду меню «Сервис — параметры», выбрать вкладку «Вид» и отметить в параметрах окна флажком пункт «формулы».
5. Теперь остается выполнить обратный ход метода Гаусса. Для этого в ячейках F17:F20 введем формулы:
| =E17-F20*D17-F19*C17-F18*B17 |
| =E18-F20*D18-F19*C18 |
| =E19-F20*D19 |
| =E20 |
В ячейках F17:F20 получим ответ:
x1 = –1,1677; x2 = 2,2446; x3 = –1,7081; x2 =2,6746.
Таблица 3.1
| A | B | C | D | E | |
| -4 | |||||
| -4 | -4 | ||||
| -8 | |||||
| =A1/$A1 | =B1/$A1 | =C1/$A1 | =D1/$A1 | =E1/$A1 | |
| =A2-A$5*$A2 | =B2-B$5*$A2 | =C2-C$5*$A2 | =D2-D$5*$A2 | =E2-E$5*$A2 | |
| =A3-A$5*$A3 | =B3-B$5*$A3 | =C3-C$5*$A3 | =D3-D$5*$A3 | =E3-E$5*$A3 | |
| =A4-A$5*$A4 | =B4-B$5*$A4 | =C4-C$5*$A4 | =D4-D$5*$A4 | =E4-E$5*$A4 | |
| =A5:E5 | =A5:E5 | =A5:E5 | =A5:E5 | =A5:E5 | |
| =B6/$B6 | =C6/$B6 | =D6/$B6 | =E6/$B6 | ||
| =B7-B$10*$B7 | =C7-C$10*$B7 | =D7-D$10*$B7 | =E7-E$10*$B7 | ||
| =B8-B$10*$B8 | =C8-C$10*$B8 | =D8-D$10*$B8 | =E8-E$10*$B8 | ||
| =A9:E10 | =A9:E10 | =A9:E10 | =A9:E10 | =A9:E10 | |
| =A9:E10 | =A9:E10 | =A9:E10 | =A9:E10 | =A9:E10 | |
| =C11/$C11 | =D11/$C11 | =E11/$C11 | |||
| =C12-C$15*$C12 | =D12-D$15*$C12 | =E12-E$15*$C12 | |||
| =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | |
| =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | |
| =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | =A13:E15 | |
| =D16/$D16 | =E16/$D16 |
В таблице 3.2 приведены результаты решения системы (после снятия флажка с пункта «формулы» меню «Сервис — параметры»).
Таблица 3.2
| A | B | C | D | E | F | |
| -4 | 12,0000 | |||||
| -4 | -4 | 14,0000 | ||||
| -8 | 20,0000 | |||||
| 10,0000 | ||||||
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 | ||
| 0,0000 | 8,0000 | -2,0000 | -5,0000 | 8,0000 | ||
| 0,0000 | 0,5000 | -5,0000 | 0,5000 | 11,0000 | ||
| 0,0000 | 0,0000 | 9,0000 | 5,0000 | -2,0000 | ||
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 | ||
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 | ||
| 0,0000 | -4,8750 | 0,8125 | 10,5000 | |||
| 0,0000 | 9,0000 | 5,0000 | -2,0000 | |||
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 | ||
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 | ||
| 1,0000 | -0,1667 | -2,1538 | ||||
| 0,0000 | 6,5000 | 17,3846 | ||||
| 1,0000 | 0,5000 | -1,0000 | 0,5000 | 3,0000 | -1,1677 | |
| 0,0000 | 1,0000 | -0,2500 | -0,6250 | 1,0000 | 2,2446 | |
| 0,0000 | 0,0000 | 1,0000 | -0,1667 | -2,1538 | -1,7081 | |
| 1,0000 | 2,6746 | 2,6746 |
6. Для проверки правильности введенных формул умножим исходную матрицу коэффициентов, хранящихся в ячейках A1:D4, на столбец найденных решений F17:F20, и запишем полученные значения в столбец F1:F4.
Для этого выделим диапазон F1:F4, введем формулу
=МУМНОЖ(A1:D4;F17:F20)
и, удерживая нажатыми клавиши Ctrl и Shift, нажмем Enter. Значения в столбцах E1:E4 и F1:F4 должны совпадать.
Замечание.Приведенный лист с формулами можно применить для решения любой системы уравнений с четырьмя неизвестными. Для этого в диапазоне A1:E4 введем коэффициенты и правые части новой системы уравнений. В диапазоне F17:F20 автоматически отобразится решение новой системы уравнений.
Проверить правильность алгоритма можно, например, заменив числа в столбце E1:E4 суммами коэффициентов уравнений. Тогда вектор решений должен состоять из единиц.
Приведем программу на языке C++ для решения системы линейных уравнений методом Гаусса с выбором главного элемента:
#include <except.h>
#include <iostream.h>
int gauss(long double **a, long double *b, long double *x, const int n);
int main(){
long double **a; long double *b; long double *x; int i,j,n;
cout <<" input n = " ; cin >> n;
try {
a= new long double*[n]; for(i=0;i<n;i++) a[i]=new long double[n];
b= new long double[n]; x= new long double[n];
}
catch (xalloc){cout <<" Could not allocate "; exit(-1);}
cout <<" input matrix a ";
for (i=0; i<n; i++)for (j=0; j<n; j++)cin >> a[i][j];
cout <<" input vector b ";
for (i=0; i<n; i++)cin >> b[i];
cout << " matrix a: ";
for (i=0; i<n; i++){cout << " "; for (j=0; j<n; j++) cout <<" "<< a[i][j];}
cout << " vector b: ";
for (i=0; i<n; i++)cout << " " << b[i];
gauss(a, b, x, n);
cout << " vector x: ";
for (i=0; i<n; i++)cout << " x[" << i <<"] = " << x[i];
cout << " Press Key and Enter to continue";
cin >> i; // for pause
for(i=0;i<n;i++)delete[] a[i];
delete a;
delete[] b;
delete[] x;
return 0;
}//end main
int gauss(long double **a, long double *b, long double *x, const int n){
// Метод Гаусса с выбором главного элемента
int i, k, m, im; long double amm, aim, r;
for (m = 0; m <= n-2; m++) {// m
amm = a[m][m]; im = m;
for (i = m; i <= n-1; i++)
if(a[i][m] > amm){amm = a[i][m]; im = i; }
if(im != m){r = b[im]; b[im] = b[m]; b[m] = r;
for (k = m; k <= n-1; k++)
{r = a[im][k]; a[im][k] = a[m][k]; a[m][k] = r;}
}
for (k = m; k <= n-1; k++)a[m][k] = a[m][k]/amm; // 3.16
b[m] = b[m] / amm; //
for (i = m + 1; i <= n-1; i++){// i
aim = a[i][m];
for (k = m; k <= n-1; k++)
a[i][k] = a[i][k] - a[m][k]*aim; // 3.17
b[i] = b[i] - b[m]*aim; //
}// end i
}// end m
x[n-1] = b[n-1]/a[n-1][n-1]; // 3.19
for (i = n - 2; i >= 0; i--){// i
x[i] = b[i]; //
for (k = i + 1; k < n; k++) // 3.20
x[i] = x[i] - a[i][k]*x[k]; //
}// end i
return 0;
}// end gauss
Приведем для сравнения вариант процедуры обычного метода Гаусса:
int gauss(long double **a, long double *b, long double *x, const int n){
// Метод Гаусса
int i, k, m; long double amm, aim;
for (m = 0; m <= n-2; m++) {// m
amm = a[m][m];
for (k = m; k <= n-1; k++)a[m][k] = a[m][k]/amm; // 3.16
b[m] = b[m] / amm;
for (i = m + 1; i <= n-1; i++){// i
aim = a[i][m];
for (k = m; k <= n-1; k++)
a[i][k] = a[i][k] - a[m][k]*aim; // 3.17
b[i] = b[i] - b[m]*aim;
} // end i
} // end m
x[n-1] = b[n-1]/a[n-1][n-1]; // 3.19
for (i = n - 2; i >= 0; i--){// i
x[i] = b[i];
for (k = i + 1; k < n; k++) // 3.20
x[i] = x[i] - a[i][k]*x[k];
} // end i
return 0;
} // end gauss
Отметим, что в процедуре gauss() массивы исходных данных a и b не сохраняются, так как результаты преобразований записываются в эти же массивы.