3.1. Метод исключения гаусса для систем линейных уравнений
Из курса линейной алгебры [Крылов и др., 1972; Курош, 1962; Фадеев, Фадеева, 1963 и др.] известно, что решение системы линейных уравнений можно просто найти по правилу Крамера- через отношение определителей. Но этот способ не очень удобен для решения систем уравнений с числом неизвестных > 5, т.е. когда найти определитель сложно, а при числе неизвестных > 10 нахождение определителя с достаточно высокой степенью точности становитсясамостоятельной вычислительной задачей. В этих случаях применяют иные методы решения, среди которых самым распространеннымявляетсяметод Гаусса.
Запишем систему линейных уравнений (1.22) в виде
(1.24)
Если матрица системы верхняя треугольная, т.е. ее элементы ниже главной диагонали равны нулю, то все хjможно найти последовательно, начиная с хn, по формуле
.
(1.25)
При j > kи аjj № 0 этот метод дает возможность найти решение системы.
Метод Гаусса для произвольной системы (1.22) основан на приведении матрицы А системы к верхней или нижней треугольной. Для этого вычитают из второго уравнения системы первое, умноженное на такое число, при которома21 = 0, т.е. коэффициент прих1во второй строке должен быть равен нулю. Затем таким же образом исключают последовательноа31 , а41 , ..., аm1 .После завершения вычислений все элементыпервого столбца, кромеа11, будут равны нулю.
Продолжая этот процесс, исключают из матрицы А все коэффициентыаij, лежащие ниже главной диагонали.
Построим общие формулы этого процесса. Пусть исключены коэффициенты изk- 1 столбца. Тогда ниже главной диагонали должны остаться уравнения с ненулевыми элементами:
Умножим k-ю строку на число
и вычтем ее из m-й строки. Первый ненулевой элемент этой строки обратится в нуль, а остальные изменятся по формулам
(1.26)
Произведя вычисления по формулам (1.26) при всех указанных индексах, исключим элементы k-го столбца. Такое исключение назовемциклом, а выполнение всех циклов назовем прямым ходомисключения.
После выполнения всех циклов получим верхнюю треугольную матрицу Асистемы (1.24), которая легко решается обратным ходом по формулам (1.25). Если при прямом ходе коэффициентаjjоказался равен нулю, то перестановкой строк перемещаем на главную диагональ ненулевой элемент и продолжаем расчет.
Если предположить, что аjj №0, то тогда можно применить подпрограммуgaus,решающую систему изnлинейных уравнений. Это одна из подпрограмм, традиционно использующаяся в библиотеках стандартных программ для ЕС-ЭВМ. Ее перевод на языкPASCALвыполнен авторами.
Формальные параметры процедуры.Входные:n(типinteger) - порядок системы;а (типreal) - массив коэффициентов системы размеромn´n;b(типreal) - массив-строка свободных членов.Выходные:х(типreal) - массив-строка, в который помещается решение системы.
procedure gaus (n:integer; a : mas; b : mas1;
var x : mas1);
type mst = array [1..n] of real;
mss = array [1..n] of mst;
var a1 : mss; b1 : mst;
i, j, m, k : integer; h : real;
begin
for i := 1 to n do
begin
b1[i] := b[i];
for j := 1 to n do
a1 [i,j] := a[i,j];
end;
for i := 1 to n-1 do
for j := i+1 to n do
begin
a1[j,i] :=- a1[j,i] / a1[i,i];
for k := i+1 to n do
a1[j,k] := a1 [j,k] + a1[j,i]*a1[i,k];
b1[j] := b1[j] + b1[i]*a1[j,i];
end;
x[n] := b1[n] / a1[n,n];
for i := n-1 downto 1 do
begin
h := b1[i];
for j := i+1 to n do
h := h - x[j]*a1[i,j];
x[i] := h / a1[i,i];
end;
end.
Если контроль за невырожденностью матрицы Анеобходим, то можно предложить воспользоваться другой процедуройgaus1.
В отличие от первой процедуры здесь в выходных параметрах есть переменная tol, возвращающая 0 при нормальном завершении работы процедуры, или 1, если на главной диагонали один из элементов равен 0, или 2, если матрицаА размерностью больше, чем 50´50.
procedure gaus1 (n:integer; a : mas; b : mas1;
var x : mas1; var tol : integer);
type mst = array [1..50] of real;
mss = array [1..50] of mst;
var a1 : mss; b1 : mst;
i, j, m, k : integer; h : real;
begin
for i := 1 to n do
begin
b1[i] := b[i];
for j := 1 to n do a1 [i,j] := a[i,j];
end;
tol := 0;
if n > 50 then
begin
tol := 2;
exit;
end;
for i := 1 to n-1 do
for j := i+1 to n do
begin
if abs(a1[i,i])< 1.0e-10 then
begin
tol := 1;
exit;
end;
a1[j,i] :=- a1[j,i] / a1[i,i];
for k := i+1 to n do
a1[j,k] := a1 [j,k] + a1[j,i]*a1[i,k];
b1[j] := b1[j] + b1[i]*a1[j,i]
end;
if abs(a1[n,n]) < 1.0e-10 then
begin
tol := 1;
exit;
end;
x[n] := b1[n] / a1[n,n];
for i := n-1 downto 1 do
begin
if abs(a1[i,i]) < 1.0e-10 then
begin
tol := 1;
exit;
end;
h := b1[i];
for j := i+1 to n do h := h - x[j]*a1[i,j];
x[i] := h / a1[i,i];
end;
end.
По поводу приведенных алгоритмов можносказать, что они не самые эффективные и пригодны для решения систем, имеющих невырожденные матрицыАиВс рангом не более 50, составленные из приблизительно одинаковых коэффициентов.Если между наибольшим и наименьшим из коэффициентов расстояние более чем 104, то эти алгоритмы применять не рекомендуется, так как погрешность вычислений будет такова, что может либо привести к вырожденной матрицеА, либо дать в результате векторХс большой вычислительной погрешностью.
Для проверки работы процедур решим систему линейных уравнений методом Гаусса с точностьюдо 0.0001:
Коэффициенты, промежуточные результаты и окончательное решение приводятся в табл. 1.14. Подобные таблицы удобно использовать для ручного счета и проверки хода решения.
Таблица 1.14
-
Коэффициенты при неизвестных
Cвободные
Cтрочные
x1
x2
x3
x4
члены
суммы
0.68
0.21
0.11
-0.08
0.55
-0.13
0.84
0.15
-0.11
0.27
0.28
0.50
0.88
0.80
-0.06
0.12
2.15
0.44
0.83
1.16
2.85
-0.01
1.44
0.61
1
0.0735
-0.1618
0.1176
3.1618
4.1912
-0.1454
-0.18319
0.15590
0.30398
0.26220
-0.51290
-0.8247
0.0729
-0.1106
-0.22398
-0.48220
1.41290
-0.88015
-0.97847
0.94530
1
-2.09060
5.6719
1.54040
6.1217
-1.47643
-0.18697
4.79139
-0.99480
0.79920
1.17230
4.1136
-0.0095
1
-3.24410
-0.54110
-2.7851
-1.60130
1.07110
-0.5302
1
-0.66890
0.3311
2.8264
-0.3337
-2.7110
-0.6689
: РЕШЕНИЕ