HOME ПРИМЕРЫ THANKS НОВИЧКАМ ДОКИ LINKS JavaScript Mail |
|
|
Turbo Pascal Examples |
Графика: |
const n=5; eps=0.00001; { all numbers less than eps are equal 0 } type matr=array[1..n,1..n] of real; var a,b,a0:matr; i,j,imx,np:integer; s0,s1:real; procedure PrintMatr(m,m1:matr;n,nz,nd:integer); var i,j:integer; begin for i:=1 to n do begin if (i=1) then write(np:2,':') else write(' '); for j:=1 to n do write(m[i,j]:nz:nd); for j:=1 to n do write(m1[i,j]:nz:nd); writeln; end; inc(np); end; 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; { отладочные присвоения a[1,1]:= 3; a[1,2]:=-1; a[1,3]:= 2; a[1,4]:= 0; a[2,1]:=-2; a[2,2]:= 1; a[2,3]:= 0; a[2,4]:= 5; a[3,1]:= 1; a[3,2]:= 4; a[3,3]:=-2; a[3,4]:= 2; a[4,1]:= 0; a[4,2]:=-2; a[4,3]:= 3; a[4,4]:=-4; a[1,1]:= 5; a[1,2]:= 7; a[1,3]:= 7; a[1,4]:= 1; a[2,1]:= 6; a[2,2]:= 6; a[2,3]:= 3; a[2,4]:= 4; a[3,1]:= 5; a[3,2]:= 1; a[3,3]:= 1; a[3,4]:= 1; a[4,1]:= 3; a[4,2]:= 3; a[4,3]:= 3; a[4,4]:= 3; } for i:=1 to n do for j:=1 to n do a0[i,j]:=a[i,j]; writeln('Starting matrix:'); np:=0; PrintMatr(a,b,n,6,1); 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])); { PrintMatr(a,b,n,6,1);} { Прямой ход } 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]); { PrintMatr(a,b,n,6,1);} end else begin writeln('Обратной матрицы не существует.'); halt; end end; {writeln('Обратный ход:');} 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; { PrintMatr(a,b,n,8,4);} end else writeln('Обратной матрицы не существует.'); MultMatr(a0,b,a); writeln('Начальная матрица, обратная к ней матрица:'); PrintMatr(a0,b,n,7,3); writeln('Проверка: должна быть единичная матрица.'); PrintMatr(a,a,n,7,3); { Выполним еще проверку насколько полученная проверочная матрица близка к единичной. Сложим отдельно суммы модулей диагональных и недиагональных элементов. По диагонали должно быть n, а не по диагонали 0 } s0:=0; s1:=0; for i:=1 to n do for j:=1 to n do if (i=j) then s1:=s1+abs(a[i,j]) else s0:=s0+abs(a[i,j]); writeln('Сумма модулей диагональных элементов: ',s1); writeln('Сумма модулей недиагональных эл-тов : ',s0); end. |
Примечание 2. При решении систем линейных уравнений методом Гаусса используют так называемый метод "выбора ведущего элемента". Необходимость его обусловлена теми же причинами, что описаны в Примечании 1 - близостью к нулю одного из диагональных элементов. Только если в вышеизложенном примере использовалось прибавление к строке "нижнележащих" строк, то выбор ведущего элемента подразумевает перестановку строк в матрице таким образом, чтобы на месте диагонального элемента была строка с максимальным по модулю значением. На практике при решении систем линейных уравнений эта операция означает лишь перестановку уравнений в системе и не влияет на решение. Однако при нахождении обратной матрицы, казалось бы, так поступать нельзя - ведь будет изменен порядок следования строк. Тогда в конце нужно вернуть первоначальный порядок и все должно быть нормально. Практика показала, что это не совсем так. Первоначально когда писалась программа, порядок строк не восстанавливался. Проверка показала, что искомая матрица найдена. Попытки привести строки матрицы в первоначальный порядок, привели к неправильному решению. А вот если порядок не восстанавливать, то решение оказывается правильным. Может кто-нибудь объяснит, почему так происходит? const n=4; eps=0.00001; { all numbers less than eps are equal 0 } type matr=array[1..n,1..n] of real; vect=array[1..n] of byte; var a,b,a0:matr; v:vect; i,j,imx,np:integer; max_v,s0,s1:real; procedure PrintMatr(m,m1:matr;n,nz,nd:integer); var i,j:integer; begin for i:=1 to n do begin if (i=1) then write(np:2,':') else write(' '); for j:=1 to n do write(m[i,j]:nz:nd); for j:=1 to n do write(m1[i,j]:nz:nd); writeln(v[i]:nz); end; inc(np); end; 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; procedure SwapRows(var a,b:matr;var v:vect;n,imx,i:integer); var j,ti:byte; tr:real; begin for j:=1 to n do begin tr:=a[imx,j]; a[imx,j]:=a[i,j]; a[i,j]:=tr; tr:=b[imx,j]; b[imx,j]:=b[i,j]; b[i,j]:=tr; end; ti:=v[imx]; v[imx]:=v[i]; v[i]:=ti; end; begin randomize; for i:=1 to n do begin for j:=1 to n do begin b[i,j]:=0; a[i,j]:=random(8); end; b[i,i]:=1; v[i]:=i; end; { a[1,1]:= 3; a[1,2]:=-1; a[1,3]:= 2; a[1,4]:= 0; a[2,1]:=-2; a[2,2]:= 1; a[2,3]:= 0; a[2,4]:= 5; a[3,1]:= 1; a[3,2]:= 4; a[3,3]:=-2; a[3,4]:= 2; a[4,1]:= 0; a[4,2]:=-2; a[4,3]:= 3; a[4,4]:=-4; a[1,1]:= 5; a[1,2]:= 7; a[1,3]:= 7; a[1,4]:= 1; a[2,1]:= 6; a[2,2]:= 6; a[2,3]:= 3; a[2,4]:= 4; a[3,1]:= 5; a[3,2]:= 1; a[3,3]:= 1; a[3,4]:= 1; a[4,1]:= 3; a[4,2]:= 3; a[4,3]:= 3; a[4,4]:= 3; } for i:=1 to n do for j:=1 to n do a0[i,j]:=a[i,j]; writeln('Starting matrix:'); np:=0; PrintMatr(a,b,n,6,1); for i:=1 to n do begin { Выбор ведущего элемента } imx:=i; max_v:=abs(a[imx,i]); for j:=i+1 to n do if (abs(a[j,i])>max_v) then begin max_v:=abs(a[j,i]); imx:=j; end; { PrintMatr(a,b,n,6,1);} if (imx<>i) then SwapRows(a,b,v,n,imx,i); { PrintMatr(a,b,n,6,1);} { Прямой ход } 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]); PrintMatr(a,b,n,6,1); end; end; writeln('Returning move:'); 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; PrintMatr(a,b,n,8,4); end else writeln('Matrix doesn''t exists.'); { return rows order } {for i:=1 to n do if (v[i]<>i) then SwapRows(a,b,v,n,v[i],i);} MultMatr(a0,b,a); writeln('Starting matrix, result:'); PrintMatr(a0,b,n,8,4); writeln('check:'); PrintMatr(a,a,n,8,5); s0:=0; s1:=0; for i:=1 to n do for j:=1 to n do if (i=j) then s1:=s1+abs(a[i,j]) else s0:=s0+abs(a[i,j]); writeln('Diagonal summ:',s1,' Not diagonal summ:',s0); end.             |
HOME |