|
Расчет обратной матрицы, Программа на Delphi |
|
|
|
Oct 2 2011, 19:25
|
Местный
  
Группа: Свой
Сообщений: 447
Регистрация: 16-11-08
Из: Украина, Донецк
Пользователь №: 41 684

|
Здравствуйте! Столкнулся по учебе с написанием программы, которая делает различные операции с матрицами (умножение, деление, вычитание, транспонирование и т.д.). Все операции сделал, кроме одной - обратная матрица. В качестве метода расчета выбрал метод расчета с помощью матрицы алгебраических дополнений. Принцип расчета - находится определитель матрицы. Затем создается транспонированная матрица алгебраических дополнений, каждый элемент которой делится на определитель матрицы. Расчет определителя - получился и работает. А вот с алгебраическими дополнения возникли трудности 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 есть функции, которые позволяют рассчитывать матрицы более проще, но мне хочется разобраться с минорами.
|
|
|
|
|
 |
Ответов
(1 - 10)
|
Oct 3 2011, 07:44
|
Местный
  
Группа: Свой
Сообщений: 447
Регистрация: 16-11-08
Из: Украина, Донецк
Пользователь №: 41 684

|
Нашел пример расчета обратной матрицы по ссылке: 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.
|
|
|
|
|
Oct 3 2011, 10:50
|
Местный
  
Группа: Свой
Сообщений: 447
Регистрация: 16-11-08
Из: Украина, Донецк
Пользователь №: 41 684

|
Подумал, что все таки лучше продолжить свое начатое решение в первом посте
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". Что я сделал неправильно?
|
|
|
|
|
Oct 3 2011, 15:10
|

Гуру
     
Группа: Свой
Сообщений: 2 399
Регистрация: 10-05-06
Из: г. Новочеркасск
Пользователь №: 16 954

|
Цитата(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 (выделите часть сообщения и нажмите крайнюю правую кнопку выше окна ввода).
|
|
|
|
|
Oct 3 2011, 18:10
|
Местный
  
Группа: Свой
Сообщений: 447
Регистрация: 16-11-08
Из: Украина, Донецк
Пользователь №: 41 684

|
Подредактировал программу. Сейчас ее код имеет следующий вид: Код 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 находится верно и вроде как циклы в отладчике правильно проходят, но программа почему-то работает не так, как нужно. Что-то в ней не так...
|
|
|
|
|
Oct 4 2011, 04:57
|
Местный
  
Группа: Свой
Сообщений: 447
Регистрация: 16-11-08
Из: Украина, Донецк
Пользователь №: 41 684

|
Ниже есть другие "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;, что он был в цикле?
|
|
|
|
|
Oct 4 2011, 15:07
|

Ambidexter
    
Группа: Свой
Сообщений: 1 589
Регистрация: 22-06-06
Из: Oxford, UK
Пользователь №: 18 282

|
1) Всё уже сделано до нас. Найдите книгу Дж.Форсайт, М.Малькольм, К.Моулер "Машинные методы математических вычислений", Мир, 1980. В параграфе 3.3 приведено описание и текст программ DECOMP и SOLVE. Вам нужна первая. Вся программа занимает 3 листа, из них половина - комментарии.
2) Определитель считать не кошерно, более показательным параметром является число обусловленности матрицы cond(A). В книге оно тоже описывается (п.3.2). Если интересуетесь, то в конце книги есть обсуждение и программа сингулярного разложения матриц, весьма мощное вычислительное средство для анализа матричных задач (п.9.2).
--------------------
Делай сразу хорошо, плохо само получится
|
|
|
|
|
  |
2 чел. читают эту тему (гостей: 2, скрытых пользователей: 0)
Пользователей: 0
|
|
|