Алгоритм метода Гаусса для обратной матрицы
- Нужно записать рядом две матрицы - слева данную матрицу, а справа единичную.
- Дальнейшие преобразования проводить как для левой матрицы, так и для правой одновременно.
- Умножить первую строку исходной матрицы на число обратное первому элементу. В результате элемент первой строки станет равным единице, а для второй матрицы - обратному первому элементу первой строки исходной матрицы.
- Умножить первую строку на число, равное первому элементу второй строки. Первый элемент первой строки станет равным первому элементу второй строки. Для второй матрицы первый элемент будет равен отношению первых элементов второй и первой строки исходной матрицы.
- Вычесть из второй строки первую строку. На первом месте второй строки окажется 0. Первый элемент второй строки правой матрицы станет равным противоположному для первого элемента первой строки.
- Вычесть из третьей и всех последующих строк первую строку, предварительно умножив ее на элемент уменьшаемой строки. Первый столбец на этом этапе будет обнулен, а в правой матрице в первом столбце ниже первой строки появятся отличные от нуля элементы.
- Вычеркнуть первую строку и первый столбец из полученных матриц на 6 этапе.
- К полученной матрице на единицу меньшего порядка применить действия по пунктам 3-7 Метода Гаусса не забывая при этом производить указанные действия над полными строками обеих матриц. В результате первый столбец правой матрицы уже будет алгебраической суммой двух слагаемых.
- Повторить пункт 9 ко всем подматрицам более низкого порядка. В результате выполнения пунктов 1-11 мы получим матрицу ниже главной диагонали которой все элементы равны нулю, тем не менее в правой матрице элементы ниже главной диагонали станут отличными от нуля.
- Выполнить аналогичные действия по обнулению элементов левой матрицы в обратном порядке снизу вверх.
- Разделить все строки матриц на значения диагональных элементов.
Слева получилась единичная матрица, а справа обратная матрица.
Решение можно разбить на два этапа: "прямой" и "обратный" ход.
Есть исходная матрица, приписываем к ней единичную:
а[1,1], а[1,2], ..., а[1,n]; 1, 0, ..., 0;
а[2,1], а[2,2], ..., а[2,n]; 0, 1, ..., 0;
...
а[n,1], а[n,2], ..., а[n,n]; 0, 0, ..., 1;
Во время "прямого" хода необходимо добиться нулей под главной диаганалью левой матрицы. Берем первую строку и делим ее на а[1,1]. Теперь на месте а[1,1] стоит 1. Вычитаем из второй строки первую умноженную на а[2,1] - на месте этого эл-та образуется ноль. Аналогично для всех строк до n-ной. Теперь в первом столбце матрицы ниже единицы стоят нули.
Переходим ко второй строке и для всех строк ниже второй повторям описанную процедуру. Теперь ниже диагонали и во второй строке нули. Так продолжаем до (n-1)-ой строки. Для n-ной строки достаточно разделить ее на а[n,n]. Матрица а приведена к верхней треугольной. На месте единичной образовалась некая матрица.
Примечание 1. Если на месте диагонального элемента левой матрицы образуется число близкое к нулю, то деление на малое число приведет к значительной погрешности вычисления. Поэтому необходимо, чтобы это число было "далеко" от нуля. С этой целью предпринимается следующий шаг: перед тем, как разделить строку на этот элемент, прибавим к строке все нижележащие строки (умноженные на -1 если в этом столбце стоит отрицательный элемент).
Обратный ход. Здесь сначала добиваемся нулей в последнем столбце матрицы а. Для этого из каждой строки (i) выше n-ной вычитаем n-ную умноженную на а[i,n]. Затем добиваемся нулей в (n-1)-ом столбце и так далее до второго столбца.
Все. Теперь слева имеем единичную матрицу, а справа, на месте единичной - искомая обратная матрица. Для проверки перемножим ее на начальную - должна получиться единичная.
const n=5;
eps=0.00001;
type matr=array[1..n,1..n] of real;
var a,b,a0:matr;
i,j,imx,np:integer;
s0,s1:real;
procedure MultString(var a,b:matr;i1:integer;r:real);
var j:integer;
begin
for j:=1 to n do
begin
a[i1,j]:=a[i1,j]*r;
b[i1,j]:=b[i1,j]*r;
end;
end;
procedure AddStrings(var a,b:matr;i1,i2:integer;r:real);
{ Процедура прибавляет к i1 строке матрицы a i2-ю умноженную на r}
var j:integer;
begin
for j:=1 to n do
begin
a[i1,j]:=a[i1,j]+r*a[i2,j];
b[i1,j]:=b[i1,j]+r*b[i2,j];
end;
end;
procedure MultMatr(a,b:matr;var c:matr);
var i,j,k:byte;
s:real;
begin
for i:=1 to n do
for j:=1 to n do
begin
s:=0;
for k:=1 to n do
s:=s+a[i,k]*b[k,j];
c[i,j]:=s;
end;
end;
function sign(r:real):shortint;
begin
if (r>=0) then sign:=1 else sign:=-1;
end;
begin { начало основной программы }
randomize; { используем автозаполнение матрицы случайными числами }
for i:=1 to n do
begin
for j:=1 to n do
begin
b[i,j]:=0;
a[i,j]:=1.0*random(8)-4;
end;
b[i,i]:=1;
end;
for i:=1 to n do
for j:=1 to n do
a0[i,j]:=a[i,j];
writeln('Starting matrix:'); np:=0;
PrintMatr( );
for i:=1 to n do
begin
{ К i-той строке прибавляем (или вычитаем) j-тую строку
взятую со знаком i-того элемента j-той строки. Таким образом,
на месте элемента a[i,i] возникает сумма модулей элементов i-того
столбца (ниже i-той строки) взятая со знаком бывшего элемента a[i,i],
равенство нулю которой говорит о несуществовании обратной матрицы }
for j:=i+1 to n do
AddStrings(a,b,i,j,sign(a[i,i])*sign(a[j,i]));
{ Прямой ход }
if (abs(a[i,i])>eps) then
begin
MultString(a,b,i,1/a[i,i]);
for j:=i+1 to n do
AddStrings(a,b,j,i,-a[j,i]);
end
else
begin
writeln('Обратной матрицы не существует.');
halt;
end
end;
if (a[n,n]>eps) then
begin
for i:=n downto 1 do
for j:=1 to i-1 do
begin
AddStrings(a,b,j,i,-a[j,i]);
end;
end
else writeln('Обратной матрицы не существует.');
MultMatr(a0,b,a);
writeln('Начальная матрица, обратная к ней матрица:');
PrintMatr( );
writeln('Проверка: должна быть единичная матрица.');
PrintMatr();
end.