Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Расчет обратной матрицы
Форум разработчиков электроники ELECTRONIX.ru > Сайт и форум > В помощь начинающему > Программирование
Andbiz
Здравствуйте! Столкнулся по учебе с написанием программы, которая делает различные операции с матрицами (умножение, деление, вычитание, транспонирование и т.д.). Все операции сделал, кроме одной - обратная матрица.

В качестве метода расчета выбрал метод расчета с помощью матрицы алгебраических дополнений.

Принцип расчета - находится определитель матрицы. Затем создается транспонированная матрица алгебраических дополнений, каждый элемент которой делится на определитель матрицы.
Расчет определителя - получился и работает. А вот с алгебраическими дополнения возникли трудности

det:=1.0;
// начинаю прямой ход Гаусса
for k:=1 to n do
begin
det:=det*A[k,k]; //вычисление определителя
for j:=k+1 to n do
begin
A[k,j]:=A[k,j]/A[k,k];
end;
for i:=k+1 to n do //начало вложенного цикла
for j:=k+1 to m do
begin
r:=A[k,j]*A[i,k];
A[i,j]:=A[i,j]-r;
end;
end;
// Делаю транспонирование матрицы
begin
for i:=1 to n do
for j:=1 to m do
D[i,j]:=A[j,i];
for i:=1 to n do
for j:=1 to m do
C[i,j]:=(1/det)*D[i,j];
begin
for i:=1 to n do
for j:=1 to m do
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу полученное значение в матрицу С
end;
end;


Красным отмечен участок, который я неправильно написал. Я просто транспонировал матрицу (заменил строки столбцами) и разделил каждый элемент на определитель, но это оказалось неправильным расчетом.

Вот пример расчета обратной матрицы.
Нажмите для просмотра прикрепленного файла

Здесь находят так называемые миноры, которые и образуют обратную матрицу.
В интернете много примеров расчета обратных матриц на паскале, но я не могу их понять и поэтому решил писать свою программу. Определитель посчитан, нужно разобраться как находить миноры, а потом их уже домножать как в примере на -1 в степени суммы индексов и делить на определитель, расположив в новой матрице.
В интернете это описывают так:
Цитата
"Создадим массив для хранения союзной матрицы. Думаю нужно создать цикл в цикле (перебор всех членов матрицы). можно создать еще один, временный массив в котором будем хранить минор. ну и все опять же вызываем функцию по нахождению определителя и вписываем результат в союзную матрицу."

Как написать код, чтобы найти эти миноры?

P.S.
Я знаю, что в Delphi есть функции, которые позволяют рассчитывать матрицы более проще, но мне хочется разобраться с минорами.
Andbiz
Нашел пример расчета обратной матрицы по ссылке:

http://www.cyberforum.ru/pascal/thread29784.html

Пытаюсь эту чужую программу вставить в свою программу.
В ней заранее предопределены функции в начале программы. У меня при попытки компилирования (поиске ошибок) пишется ошибка перед функциями Det и AlgDop в теле программы:
Цитата
43. Statement expected but Function found. В этом месте должен стоять оператор. Сообщение выдается во всех случаях, когда в тело блока или секцию инициализации ошибочно помещают описание (<Что-то>). Ошибочная форма обращения к процедуре Procedure <Имя> или к функции Function <Имя> также вызывает сообщение.


У меня расчет матрицы происходит по "многоступенчатому" условию if, then, в каждую из ступенек которого заложен определенный расчет матрицы. Если функцию убрать (function Det), то возле некоторых переменных внутри функций появляется ошибка

Цитата
36. Not enough actual parameters. He хватает фактических параметров.


Т.е. в моем понимании программа ранее обращалась к переменным, которые были расчитаны после формулы (предварительная функция), то сейчас расчет идет сверху вниз и программе не хватает данных для расчета. А объявить функцию у меня не получается, т.к. она должна находиться в цикле, а Делфи это почему-то не нравится. Как быть в этом случае?

Текст программы
Цитата
procedure TForm1.Button_ClearClick(Sender: TObject);
var i,j:integer;
begin
Edit_n.Clear; // Очищаю значение Edit_n
Edit_m.Clear; // Очищаю значение Edit_m
with StringGrid_A do // Очищаю значения элементов матрицы А
for i:=0 to RowCount do
for j:=0 to ColCount do
Cells[j,i]:='';
with StringGrid_B do // Очищаю значения элементов матрицы B
for i:=0 to RowCount do
for j:=0 to ColCount do
Cells[j,i]:='';
with StringGrid_C do // Очищаю значения элементов матрицы C
for i:=0 to RowCount do
for j:=0 to ColCount do
Cells[j,i]:='';
end;

procedure TForm1.Button_CloseClick(Sender: TObject);
begin
close
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
Randomize // Генерация случайных чисел при запуске программы
end;

procedure TForm1.Button_CalcClick(Sender: TObject);
const n_max=15; // максимальное число строк матриц равно 15
const m_max=15; // максимальное число столбцов матриц равно 15
type matr_1=array [1..n_max, 1..m_max] of real; // числа строк и столбцов матриц - любые числа
type matr_2=array [1..n_max, 1..m_max] of real; // числа строк и столбцов матриц - любые числа
function Det(var size:integer;e:matr_2):real;forward;
function AlgDop(var size:integer;e:matr_2):real;forward;
var n,m:integer; // число строк и столбцов - целое число
var a,e,b:matr_1; // a, b - это матрицы, элеметы которых любые числа
var temp,c,d:matr_2; // c - это матрица с возможными дробями в качестве элементов матрицы
size,i,j,k:integer; // индексы i и j - целые числа
r,s,dt:real; // det - определитель системы, r - произведение при прямом ходе

begin
// Ввод размерности матриц и настройка компонентов StringGrid (матриц)
n:=StrToInt(Edit_n.Text); // Принимаю число строк равным n
m:=StrToInt(Edit_m.Text); // Принимаю число строк равным m
StringGrid_A.RowCount:=n; // для матрицы А число строк принимаю равным n
StringGrid_A.ColCount:=m; // для матрицы А число столбцов принимаю равным m
StringGrid_B.RowCount:=n; // для матрицы B число строк принимаю равным n
StringGrid_B.ColCount:=m; // для матрицы B число столбцов принимаю равным m
StringGrid_C.RowCount:=n; // для матрицы C число строк принимаю равным n
StringGrid_C.ColCount:=m; // для матрицы C число столбцов принимаю равным m

// Ввод матрицы А
if (RadioGroup1.ItemIndex=0) // заполнение вручную
then for i:=1 to n do
for j:=1 to m do
A[i,j]:=StrToInt(StringGrid_A.Cells[j-1,i-1])
else for i:=1 to n do // заполнение случайными числами
for j:=1 to m do
begin
A[i,j]:=Random(101); // случайное число от 0 до 100
StringGrid_A.Cells[j-1,i-1]:=FloatToStr(A[i,j]) // вношу полученное значение в матрицу А
end;

// Ввод матрицы B
if (RadioGroup2.ItemIndex=0) // заполнение вручную
then for i:=1 to n do
for j:=1 to m do
B[i,j]:=StrToInt(StringGrid_B.Cells[j-1,i-1])
else for i:=1 to n do // заполнение случайными числами
for j:=1 to m do
begin
B[i,j]:=Random(101); // случайное число от 0 до 100
StringGrid_B.Cells[j-1,i-1]:=FloatToStr(B[i,j]) // вношу полученное значение в матрицу B
end;

// Заполнение матрицы С
if (RadioGroup3.ItemIndex=0) // сложить матрицы (А+В)
then for i:=1 to n do
for j:=1 to m do
begin
C[i,j]:=A[i,j]+B[i,j]; // процедура сложения элементов матриц
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу полученное значение в матрицу С
end
else
if (RadioGroup3.ItemIndex=1) // вычесть матрицы (А-В)
then for i:=1 to n do
for j:=1 to m do
begin
C[i,j]:=A[i,j]-B[i,j]; // процедура вычитания элементов матриц
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу полученное значение в матрицу С
end
else
if (RadioGroup3.ItemIndex=2) // перемножить матрицы (А*В)
then for i:=1 to n do
for j:=1 to m do
begin // начинаю цикл расчета элемента С [i,j]
s:=0;
for k:=1 to n do
s:=s+A[i,k]*B[k,j];
C[i,j]:=s;
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу полученное значение в матрицу С
end
else
if (RadioGroup3.ItemIndex=3) // транспонировать матрицу А
then for i:=1 to n do
for j:=1 to m do
begin
StringGrid_C.Cells[j-1,i-1]:=StringGrid_A.Cells[i-1,j-1]; // замена элемента строки на противоположный элемент столбцами
end
else
if (RadioGroup3.ItemIndex=4) // транспонировать матрицу B
then for i:=1 to n do
for j:=1 to m do
begin
StringGrid_C.Cells[j-1,i-1]:=StringGrid_B.Cells[i-1,j-1]; // замена элемента строки на противоположный элемент столбцами
end
else
if (RadioGroup3.ItemIndex=5) // вычисляю обратную матрицу А
then

begin
for i:=1 to n-1 do
begin
for j:=1 to m-1 do
temp[i,j]:=e[i,j];
for j:=m+1 to size do
temp[i,j-1]:=e[i,j];
end;
for i:=n+1 to size do
begin
for j:=1 to m-1 do
temp[i-1,j]:=e[i,j];
for j:=m+1 to size do
temp[i-1,j-1]:=e[i,j];
end;
if (n+m) mod 2 = 0
then AlgDop:=det(size-1,temp)
else AlgDop := -det(size-1,temp);
end;



procedure ObrMatr(a:matr;var am:matr;n:integer);
var i,j:integer;
dt:real;
begin
dt:=det(n,a);
if abs(dt)<1e-15 then begin
writeln('Malenkiy opredelitel!');
exit;
end;
for i:=1 to n do begin
for j:=1 to n do am[i,j]:=AlgDop(n,a,j,i)/dt;
end;
end;

procedure print_matr(n:integer;var a:matr);
var i,j:integer;
begin
for i:=1 to n do
begin
for j:=1 to n do
write(a[i,j]:7:3);
writeln;
end;
end;

var a,obr:matr;
i,j,n:integer;

begin
clrscr;
write('Vvedite N: ');
readln(n);
writeln('‚ўҐ¤ЁвҐ ',n*n,'elementov matricy:');
for i:=1 to n do
for j:=1 to n do
begin
write('a[',i,',',j,']=');
readln(a[i,j]);
end;
clrscr;
writeln('Ishodnaya matrica:');
print_matr(n,a);
ObrMatr(a,obr,n);
writeln('Obratnaya matrica:');
print_matr(n,obr);
readln
end.
end;
end;
end;
end.
Andbiz
Подумал, что все таки лучше продолжить свое начатое решение в первом посте

begin
det:=1.0;
// начинаю прямой ход Гаусса
for k:=1 to n do
begin
det:=det*A[k,k];//вычисление определителя
for j:=k+1 to n do
begin
A[k,j]:=A[k,j]/A[k,k];
end;
for i:=k+1 to n do//начало вложенного цикла
for j:=k+1 to m do
begin
r:=A[k,j]*A[i,k];
A[i,j]:=A[i,j]-r;
end;
end;
Edit1.Text:=FloatToStr(det);
// Делаю транспонирование матрицы А
begin
for i:=1 to n do // начинаю цикл обхода всех строк и столбиков матрицы А
for j:=1 to m do
begin
for i_1:=1 to n do // начинаю второй цикл обхода матрицы А, в котором будут отсеиваться одна строка и один столбик
for j_1:=1 to m do
if (i_1<>i) and (j_1<>j)
then
for l:=1 to n-1 do // заполняю матрицу D, считая, что количество строк алгебраического дополнения на одну меньше, чем матрицы А
for p:=1 to m-1 do // заполняю матрицу D, считая, что количество столбцов алгебраического дополнения на одного меньше, чем матрицы А
A[i,j]:=D[l,p]; // вношу соответствующий элемент в временную матрицу D
end;
begin
det_1:=1.0;
// начинаю прямой ход Гаусса
for k:=1 to n-1 do
begin
det_1:=det_1*D[k,k]; //вычисление определителя
for j:=k+1 to n-1 do
begin
D[k,j]:=D[k,j]/D[k,k];
end;
for i:=k+1 to n-1 do //начало вложенного цикла
for j:=k+1 to m-1 do
begin
r:=D[k,j]*D[i,k];
D[i,j]:=D[i,j]-r;
end;
end;
C[i,j]:=((i+j)*ln(-1)*det_1)/det; // рассчитываю соответствующий элемент матрицы С
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу рассчитанный элемент в матрицу
end;
end;
end;

При запуске в тексте, выделенном красным цветом появляется ошибка "Invalid floating point operation". Что я сделал неправильно?
Палыч
Цитата(Andbiz @ Oct 3 2011, 14:50) *
Что я сделал неправильно?


ln(-1) - логорифм натуральный от минус единицы? laughing.gif Из-за этого ошибка при выполнении.
Andbiz
Да, точно. Я тут ошибся - логарифма от -1 нет. Возвел в степень через Power. Теперь при попытке вычисления обратной матрицы 3*3 в результирующей матрице заполняется числом только элемент с индексом [3,3], все остальные без изменения. Да и то число, которое ставится - неверное (к примеру нужно 1, а программа ставит -3).
Палыч
Цитата(Andbiz @ Oct 3 2011, 16:26) *
...в результирующей матрице заполняется числом только элемент с индексом [3,3], все остальные без изменения.

Это - естественно. Последние две строчки программы
Код
C[i,j]:=((i+j)*ln(-1)*det_1)/det; // рассчитываю соответствующий элемент матрицы С
StringGrid_C.Cells[j-1,i-1]:=FloatToStr(C[i,j]) // вношу рассчитанный элемент в матрицу

расположены вне всяких циклов... Чтобы ясно была видна структура программы, обычно, используют отступы (например, соответствующий begin'у end начинают в одной и той же позиции строки, а начало строк между ними сдвигают вправо на несколько позиций). И Вы используйте отступы, а тест программы в посте обрамляйте тегами code (выделите часть сообщения и нажмите крайнюю правую кнопку выше окна ввода).
Andbiz
Подредактировал программу. Сейчас ее код имеет следующий вид:

Код
              
                         begin
                         det:=1;
                         // начинаю прямой ход Гаусса
                         for k:=1 to n do
                         begin
                              det:=det*A[k,k]; //вычисление определителя
                              for j:=k+1 to n do
                              begin
                                   A[k,j]:=A[k,j]/A[k,k];
                              end;
                              for i:=k+1 to n do //начало вложенного цикла
                              for j:=k+1 to m do
                              begin
                              r:=A[k,j]*A[i,k];
                              A[i,j]:=A[i,j]-r;
                              end;
                         end;
                         Edit1.Text:=FloatToStr(det);
                         // Создаю союзную матрицу D
                         begin
                              for i:=1 to n do  // начинаю цикл обхода всех строк и столбиков матрицы А
                              for j:=1 to m do
                              begin
                                   l:=1;
                                   p:=1;
                                   for i_1:=1 to n do // начинаю второй цикл обхода матрицы А, в котором будут отсеиваться одна строчка и один столбик
                                   for j_1:=1 to m do
                                       if (i_1<>i) and (j_1<>j)  // отсеиваю элементы, которые находятся в столбике или строчке рассматриваемого элемента матрицы
                                       then
                                           begin
                                           D[l,p]:=A[i_1,j_1]; // вношу соответствующий элемент в временную матрицу D
                                           p:=p+1;  // прибавляю 1 к индексу столбика матрицы D
                                           if(p>(m-1)) // проверяю индекс столбика на крайнее значение, если оно больше (n-1) исходной таблицы, то перехожу на новую строку
                                           then
                                           begin
                                           p:=1;  // при начале новой строки столбик приравниваю единице
                                           l:=l+1;  // и на одну строчку опускаюсь вниз, прибавив к l единицу
                                           end;
                                           end;
                                               begin
                                               if (n-1)=1 then det_1:=D[1,1]  // если матрица еденичная, то определитель определить легком
                                               else
                                               if (n-1)=2 then det_1:=(D[1,1]*D[2,2]-D[1,2]*D[2,1])  // если матрица имеет размерность 2*2, то определитель можно найти по следующей формуле
                                               else
                                                    begin
                                                    det_1:=1;
                                                    // начинаю прямой ход Гаусса
                                                    for k:=1 to n-1 do
                                                        begin
                                                        det_1:=det_1*D[k,k];  //вычисление определителя
                                                        for p:=k+1 to n-1 do
                                                            begin
                                                                 D[k,p]:=D[k,p]/D[k,k];
                                                            end;
                                                            for l:=k+1 to n-1 do  //начало вложенного цикла
                                                            for p:=k+1 to m-1 do
                                                            begin
                                                                 r:=D[k,p]*D[l,k];
                                                                 D[l,p]:=D[l,p]-r;
                                                            end;
                                                        end;
                                                    end;
                                                        C[i,j]:=(power((-1),(i+j))*det_1)/det; // рассчитываю соответствующий элемент матрицы С
                                                        StringGrid_C.Cells[i-1,j-1]:=FloatToStr(C[i,j])  // вношу рассчитанный элемент в матрицу
                                               end;


Но расчет все равно не правильный. По поводу строчки C[i,j]:=((i+j)*ln(-1)*det_1)/det; находящейся вне цикла - не вижу, как она находится вне цикла... Проверяю через отладчик - строчка каждый раз после цикла присваивает элементам матрицы новые значения. Единственный их недостаток - значения неправильные.

Нажмите для просмотра прикрепленного файла

Определитель исходной матрицы Det находится верно и вроде как циклы в отладчике правильно проходят, но программа почему-то работает не так, как нужно. Что-то в ней не так...
Палыч
Цитата(Andbiz @ Oct 3 2011, 22:10) *
По поводу строчки C[i,j]:=((i+j)*ln(-1)*det_1)/det; находящейся вне цикла - не вижу, как она находится вне цикла... Проверяю через отладчик - строчка каждый раз после цикла присваивает элементам матрицы новые значения.

В таком случае: укажите те циклы по i и j, в которых по Вашему мнению находится этот оператор. То, что Ваше утверждение неверно легко проверить: обратите внимание, что ниже указанных мной строк - только один-единственный end, который соответствует концу программы.
Andbiz
Ниже есть другие "end" - я их не скопировал, т.к. полностью вся программа больше и это участок этой программы. При копировании побоялся скоприовать лишние "end"ы (этот участок в самом конце) Примерно расставил соответстствие между begin и end цифрами. По поводу цикла по i и j, к которому относится C[i,j]:=((i+j)*ln(-1)*det_1)/det; - , то это

Код
for i:=1 to n do  // начинаю цикл обхода всех строк и столбиков матрицы А
                              for j:=1 to m do


Вот программа с "end" в конце программы

Код
    
                         begin
                         det:=1;
                         // начинаю прямой ход Гаусса
                         for k:=1 to n do
                         begin
                              det:=det*A[k,k]; //вычисление определителя
                              for j:=k+1 to n do
                              begin
                                   A[k,j]:=A[k,j]/A[k,k];
                              end;
                              for i:=k+1 to n do //начало вложенного цикла
                              for j:=k+1 to m do
                              begin
                              r:=A[k,j]*A[i,k];
                              A[i,j]:=A[i,j]-r;
                              end;
                         end;
                         Edit1.Text:=FloatToStr(det);
                         // Создаю союзную матрицу D
             (9)            begin
                              for i:=1 to n do  // начинаю цикл обхода всех строк и столбиков матрицы А
                              for j:=1 to m do
                 (8)             begin
                                   l:=1;
                                   p:=1;
                                   for i_1:=1 to n do // начинаю второй цикл обхода матрицы А, в котором будут отсеиваться одна строчка и один столбик
                                   for j_1:=1 to m do
                                       if (i_1<>i) and (j_1<>j)  // отсеиваю элементы, которые находятся в столбике или строчке рассматриваемого элемента матрицы
                                       then
                              (6)             begin
                                           D[l,p]:=A[i_1,j_1]; // вношу соответствующий элемент в временную матрицу D
                                           p:=p+1;  // прибавляю 1 к индексу столбика матрицы D
                                           if(p>(m-1)) // проверяю индекс столбика на крайнее значение, если оно больше (n-1) исходной таблицы, то перехожу на новую строку
                                           then
                             (7)              begin
                                           p:=1;  // при начале новой строки столбик приравниваю единице
                                           l:=l+1;  // и на одну строчку опускаюсь вниз, прибавив к l единицу
                           (7)                end;
                         (6)                  end;
                                (5)               begin
                                               if (n-1)=1 then det_1:=D[1,1]  // если матрица еденичная, то определитель определить легком
                                               else
                                               if (n-1)=2 then det_1:=(D[1,1]*D[2,2]-D[1,2]*D[2,1])  // если матрица имеет размерность 2*2, то определитель можно найти по следующей формуле
                                               else
                                           (4)         begin
                                                    det_1:=1;
                                                    // начинаю прямой ход Гаусса
                                                    for k:=1 to n-1 do
                                        (3)                begin
                                                        det_1:=det_1*D[k,k];  //вычисление определителя
                                                        for p:=k+1 to n-1 do
                                        (2)                    begin
                                                                 D[k,p]:=D[k,p]/D[k,k];
                                      (2)                      end;
                                                            for l:=k+1 to n-1 do  //начало вложенного цикла
                                                            for p:=k+1 to m-1 do
                                      (1)                    begin
                                                                 r:=D[k,p]*D[l,k];
                                                                 D[l,p]:=D[l,p]-r;
                                     (1)                       end;
                                         (3)               end;
                                       (4)             end;
                                                        C[i,j]:=(power((-1),(i+j))*det_1)/det; // рассчитываю соответствующий элемент матрицы С
                                                        StringGrid_C.Cells[i-1,j-1]:=FloatToStr(C[i,j])  // вношу рассчитанный элемент в матрицу
                              (5)                 end;
                                   (8)   end;

                     (9)     end;
                          end;


А куда по Вашему мнению нужно скопировать C[i,j]:=((i+j)*ln(-1)*det_1)/det;, что он был в цикле?
=GM=
1) Всё уже сделано до нас. Найдите книгу Дж.Форсайт, М.Малькольм, К.Моулер "Машинные методы математических вычислений", Мир, 1980. В параграфе 3.3 приведено описание и текст программ DECOMP и SOLVE. Вам нужна первая. Вся программа занимает 3 листа, из них половина - комментарии.

2) Определитель считать не кошерно, более показательным параметром является число обусловленности матрицы cond(A). В книге оно тоже описывается (п.3.2). Если интересуетесь, то в конце книги есть обсуждение и программа сингулярного разложения матриц, весьма мощное вычислительное средство для анализа матричных задач (п.9.2).
Andbiz
Большое спасибо за ответы, но все таки попытался свое довести до конца. Ошибку нашел.
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2025 Invision Power Services, Inc.