program Matrix_Process;
(******************************************************
* Демонстрация работы с двумерными массивами *
* (матрицами) *
******************************************************)
{ uses power1, gamma1, integral;}
const
NMAX = 20;
MMAX = 20;
type
Matrix = array [1..NMAX, 1..MMAX] of extended;
var
Funk, Matr: Matrix;
n,m: integer;
Str:string;
file1,file2:string;
s1,s2:char;
(***************************************************
Процедура ввода матрицы из файла
****************************************************)
procedure EnterMatrixFromFile;
var
i,j: integer;
f:text;
begin
assign(f,'matrixdat.txt');
reset(f);
i := 1;
while not eof(f) do
begin
j := 1 ;
while not eoln(f) do
begin
read(f,Matr[i, j]);
{writeln('Matr[',i,',',j,']',Matr[i,j]);}
inc(j);
end;
readln(f); inc(i);
end;
close(f);
n:=i-1; {Є®«ЁзҐбўв® бва®Є}
m:=j-1; {Є®«ЁзҐбвў® бв®«Ўж®ў}
writeln('ђ §¬Ґа Ёб室®© ¬ ваЁжл: ‘ва®Є n=',n,',','‘в®«Ўж®ў m=',m);
end;
(******************************************************
Функция возведения действительного числа а в степень b
******************************************************)
Function Power (a, b : extended): extended; { a^b }
begin
if (a > 0) and (b<>0) then Power := exp(b * ln (a));
if (a < 0) and (b<>0) then Power := exp(b * ln (-a));
if b = 0 then Power := 1;
if a = 0 then Power := 0;
end;
(******************************************************
Процедура печати матрицы реальной размерности (n,m) на
экране, располагающая одну строку матрицы на одной строке
экрана
******************************************************)
procedure PrintMatrix (Var Matrix_: Matrix; nn, mm: integer);
var
i, j: integer;
begin
for i:=1 to nn do begin
for j:=1 to mm do
writeln('Matrix_[',i,',',j,']',Matrix_[i,j]);
writeln;
end;
end;
(******************************************************
Процедура печати матрицы реальной размерности (n,m) в файл,
располагающая одну строку матрицы на одной строке
******************************************************)
procedure PrintMatrixToFile (Var Matrix_: Matrix; Var String_: string; nn: integer);
var
i, j: integer;
g:text;
begin
assign(g,String_);
rewrite(g);
for i:=1 to nn do
begin
for j:=1 to nn do
begin
write (g,matrix_[i,j]);
write (g,' ');
end;
writeln(g);
end;
close(g);
end;
(*******************************************************
Гамма-Функция Эйлера
*******************************************************)
(******************************************************)
function Gamma(x: extended): extended;
const
pi = 3.1415926;
var
i,j : integer;
y,gam : extended;
begin { gamma function }
If x=0 then gamma:=1;
if x>0.0 then
begin
y:=x+2.0;
gam:=sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
Gamma:=gam/(x*(x+1))
end
else { x<0 }
begin
j:=0;
y:=x;
repeat
j:=j+1;
y:=y+1.0
until y>0.0;
gam:=gamma(y); { recursive call }
for i:=0 to j-1 do
gam:=gam/(x+i);
Gamma:=gam
end { x<0 }
end; { gamma function }
(******************************************************)
(**********************************************************************
Программа вычисления определенного интеграла
tk
I= S power(t,-1+m/2)*exp(-t/2) dt
tn
по методу Симпсона.
Значение определенного интеграла по
формуле Симпсона вычисляется по формуле:
ISimps=2*h/3*(0.5*Q(A)+2*Q(A+h)+Q(A+2*h)+2*Q(A+3*h)+...
+2*Q(B-h)+0.5*Q(B))
где A=0 и B - нижняя и верхняя границы интервала интегрирования,
N - число разбиений интервала интегрирования,
h=(B-A)/N, причем N должно быть четным.
***********************************************************************)
Function Simps(b:extended;m:integer):extended;
var
sum, h: extended;
N,j:Integer;
begin
N:=1000;
{if Odd(N) then N:=N+1;}
h:=(b-a)/N;
sum:=0.5*(power(b,-1+m/2)*exp(-b/2));
for j:=1 to N-1 do
sum:=sum+(j mod 2+1)*power(j*h,-1+m/2)*exp(-j*h/2);
Simps:=2*h*sum/3
end;
(******************************************************
Процедура попарной обработки строк прямоугольной матрицы
******************************************************)
procedure Funky(Var Matrix_: Matrix; nn, mm: integer);
var
i,j,k,p,StSv: integer;
Sum1,Sum2,Sum3: extended;
dd:text;
begin
assign (dd,'matres.txt');
rewrite(dd);
for k:=1 to nn do
begin
for p:=1 to nn do
begin
Sum1:=0;
Sum2:=0;
Sum3:=0;
for j:=1 to mm do
begin
StSv:=0;
Sum1:=Sum1+Matrix_[k,j];
Sum2:=Sum2+Matrix_[p,j];
end;
{writeln('Sum1=',Sum1);
writeln('Sum2=',Sum2);
readln;}
for j:=1 to mm do
begin
If ((Matrix_[k,j]+Matrix_[p,j])<>0) then
begin
If (Sum1=0) and (Sum2<>0) then
Sum3:=Sum3+((1/(Matrix_[k,j]+Matrix_[p,j]))*
(-Matrix_[p,j]/Sum2)*(-Matrix_[p,j]/Sum2));
If (Sum2=0) and (Sum1<>0) then
Sum3:=Sum3+( (1/(Matrix_[k,j]+Matrix_[p,j]))*
(Matrix_[k,j]/Sum1)*(Matrix_[k,j]/Sum1));
If Sum1*Sum2<>0 then
Sum3:=Sum3+((1/(Matrix_[k,j]+Matrix_[p,j]))*
(Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2)*
(Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2))
end
else
inc(StSv);
Sum3:=Sum3;
end;
{Funk[k,p]:=Sum1*Sum2*Sum3;}
StSv:=mm-StSv;
Funk[k,p]:= Simps(Sum1*Sum2*Sum3,StSv)/
(Power(2,StSv/2)*Gamma(StSv/2));
write(dd,Funk[k,p]); write(dd,' ');
end;
writeln(dd);
end;
close(dd);
end;
(*******************************************************
Процедура замены в файле символа точки на символ запятой
*******************************************************)
procedure ExchSymbolInFile (var String1,String2:string; var symbol_1,symbol_2:char);
var
char_:char;
g,f:text;
begin
assign (g,String1);
reset(g);
assign (f,String2);
rewrite(f);
while not eof(g) do
begin
while not eoln(g) do
begin
read(g,char_);
if char_=symbol_1 then
write(f,symbol_2)
else
write (f,char_);
end;
readln(g);
writeln(f);
end;
close(f);
close(g);
end;
begin
file1:='matrix.txt';
file2:='matrixdat.txt';
s1:='.'; s2:=',';
ExchSymbolInFile(file1,file2,s2,s1);
EnterMatrixFromFile;
Funky (Matr,n,m);
file1:='mat_res.txt';
file2:='mat_res_v.txt';
ExchSymbolInFile(file1,file2,s1,s2);
end.