Обратная матрица методом гаусса паскаль
Найти обратную матрицу (Pascal)
Вообщем задача состоит в том, что надо найти обратную матрицу, по следующему алгоритму:
1. Разбиваем матрицу на 2 треугольные
2. Находим обратные матрицы к треугольным и перемножаем их (должны получить обратную матрицу к той что нам дана)
3. Проверка.
Вообщем загвоздка сейчас у меня в следующем: для того чтобы разбить матрицу на треугольные, нада чтобы определитель матрицы не был равен нулю. Так вот, подскажите по какому алгоритму можно найти определитель матрицы (я хотел с алгебраическими дополнениями, но там пальцы сломать можно, и мне сказали что есть другой способ).
13 ответов
Не помню что такое «с алгебраическими дополнениями», но припоминается что можно привести ее к треугольному виду (ето не то что разбить на 2 треугольние, а так чтоб все что под диагоналю было равно 0), а потом перемножыть елементы диагонали.
ЗЫ. ИМХО. А прежний аватар красивее был.
когда я тестирую, вроде все нормально, но когда матрица 2х2 (порядок равен 2), то определитель что-то не считается:confused: не могу найти ошибку, помогите пожалуйста.
(Может не совсем будет роботать потому что я прямо тут и набрал с копипейстом.)
Твоя ошибка в том что после первого ифа надо было остальное в елсе взять. Но ето машинально.
А теперь концептуально. Видиш, у тебя вычисление для н=2 делается два раза (в ифе и после цыкла). Старайся избегать таких ситуаций.
Извени, но помойму я тебе не так подсказал сам метод. Я тут на бумажке пробую и не получается у меня так посчитать. Пошел учить матчасть.
PS Ето у меня руки кривые на бумажке считать, таки правильно сказал. Ложная тривога.
Завтра покажу как без вспомагательной mat11 обойтись. Сегодня не успеваю.
type
Tmatr = array[1..NMAX, 1..NMAX] of real;
const
mat: Tmatr = ((8, 4, 9, 2),
(3, 7, 4, 5),
(3, 6, 9, 8),
(1, 1, 4, 1));
var
mat1: Tmatr;
n, i, j, k, current: integer;
det, buf: real;
BEGIN
mat1 := mat;
n := 4;
det:=1;
current := 1;
while(n > current) do
begin
if mat1[current, current]=0 then
begin
k:=0;
for i:=current to n do
if mat1[i,current] <> 0 then k:=i;
if k=0 then
begin
det:=0;
break;
end;
for j:=current to n do
begin
buf:=mat1[current,j];
mat1[current,j]:=mat1[k,j];
mat1[k,j]:=buf;
end;
det:=(-1)*det;
end;
for i:=current+1 to n do
begin
buf := mat1[i, current] / mat1[current, current];
for j:=current to n do
mat1[i,j]:=mat1[i,j]-buf*mat1[current,j];
end;
current := current + 1;
end;
for i := 1 to n do
det := det * mat1[i,i];
for i := 1 to n do
begin
for j := 1 to n do
write(mat1[i,j]:10:3, ‘, ‘);
writeln;
end;
writeln(‘det=’,det:8:3);
end.
Диагональ не сводится к единичной, а остается как есть. При обработке матрицы меняется только знак дискриминаннта или он устанавливается в 0. Потом в конце бежым по диагонале.
Замечание тебе
1.
Условие лутше ставить перед обменом местами строк. Кроме того при к=0 нужно прекращать приведение матрицы, а тто потом сразу деление на 0 вылезет.
2. В паскале можно переприсвоить масив полностю, а не поелементно. Но для етого надо чтоб они были одного типа. Тоесть примерно так.
type
Tmatr = array[1..NMAX, 1..NMAX] of real;
const
mat: Tmatr = ((8, 4, 9, 2),
(3, 7, 4, 5),
(3, 6, 9, 8),
(1, 1, 4, 1));
3. Не обязательно роботать через елемент [1,1]. Сделaл через current. Так у тебя остается неиспорченая матрица. (Выводится в конце).
4. Можно убрать вспомагательную матрицу и делать все на одной.
5. Убрал проверку на розмер матрицы ровный 2. Ето конечно лишняя итерацыя цыкла, но мое мнение что надо стараться делать код более универсальным. Тоесть чтоб он обрабатывал любые входные данные одинаково. Так снижается вероятность ошибок. В етом ты уже сам убедился. 🙂
6. Основное замечание. Индексируй масивы с нуля. Ето секономит тебе кучу времени и нервов когда начнеш писать на Си подобных языках.
Надеюсь мой код будет тебе полезным.
procedure output_matrix(n:integer;mat:mat1);<Вывод матрицы>
var i,j:integer;
begin
for i:=1 to n do
begin
for j:=1 to n do
begin
write(mat[i,j]:5:5);
write(‘ ‘);
end;
writeln;
end;
writeln;
end;
procedure bermuda(var nt,vt,d:mat1; n:integer);<Разбиение матрицы на 2 треугольные>
var i,j,k:integer;
s:real;
begin
for i:=1 to n do
begin
nt[i,1]:=d[i,1];
vt[i,1]:=0;
end;
for j:=1 to n do
begin
vt[1,j]:=d[1,j]/nt[1,1];
if j>1 then nt[1,j]:=0;
end;
for i:=2 to n do
for j:=2 to n do
begin
if i>=j then
begin
s:=0;
for k:=1 to j-1 do
s:=s+nt[i,k]*vt[k,j];
nt[i,j]:=d[i,j]-s;
vt[i,j]:=0;
if i=j then vt[i,j]:=1;
end
else
begin
s:=0;
for k:=1 to i-1 do
s:=s+nt[i,k]*vt[k,j];
vt[i,j]:=(d[i,j]-s)/nt[i,i];
nt[i,j]:=0;
end;
end;
end;
procedure ont(n:integer; var mat,obnt:mat1);<нахождение обратной нижней треугольной матрицы>
var i,j,a:integer;
s:real;
begin
for i:=1 to n do
begin
for j:=1 to n do
begin
if i=j then obnt[i,j]:=1/mat[i,j];
if i j then
begin
s:=0;
for a:=j to i-1 do
begin
s:=s+mat[i,a]*obnt[a,j];
obnt[i,j]:=-1*s/mat[i,i];
end;
end;
end;
end;
end;
procedure ovt(n:integer; var mat,obvt:mat1);<нахождение обратной верхней треугольной матрицы>
var i,j,a:integer;
s:real;
begin
for i:=n downto 1 do
begin
for j:=1 to n do
begin
if i>=j then obvt[i,j]:=mat[i,j];
if j=i+1 then obvt[i,j]:=-mat[i,j];
if j>i+1 then
begin
s:=0;
for a:=i+1 to j-1 do
begin
s:=s+mat[i,a]*obvt[a,j];
obvt[i,j]:=-1*(mat[i,j]+s);
end;
end;
end;
end;
end;
procedure multi_matrix(n:integer; m1,m2:mat1; var mm:mat1);<произведение матриц>
var i,j,a:integer;
x:real;
begin
for i:=1 to n do
for j:=1 to n do
begin
x:=0;
for a:=1 to n do x:=x+m1[i,a]*m2[a,j];
mm[i,j]:=x;
end;
end;
end.
Конечно код не маленький, но кому не трудно просмотрите пожалуйста. При тестировании он для всех матриц работает нормально, кроме тех в которых самый первый элемент (т.е. с индексом [1,1]) равен нулю, для них не считает:(
Пытаюсь найти сам, но что-то не получается (это как с русским языком, свои ошибки сразу не находишь, надо отвлечся, а на сежий глаз легче увидеть, особенно если она какая-нить глупая)
