Правила форума
Прочтите прежде чем задавать вопрос!
Спойлер (Показать/Скрыть)
{$A+,B-,D+,E-,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S-,T-,V+,X+,Y+}
Program _Approc_from_file;{Аппроксимация функций методом наименьших квадратов}
{$E-,N+,G+,D+}
Uses Graph, CRT;
type
Sto = array [ 1..70 ] of Extended;
const
n = 1368 div 2 - 1;
x0 : Extended = 1.0;
xn : Extended = 4.0;
Yb : Integer = 464;
var
m, nn, k, p, s : Byte;
ndot : Word;
Dr, R, e : Integer;
y : array [0..n] of Extended;
xxx : array [0..n] of Extended;
dx, ymax, ymin, x : Extended;
A : array [ 1..70 ] of Sto;{Массив }
B : Sto ;
i, j : Longint;
ff : Text;
Procedure Swap(p, q: Byte);
var
j : Byte;
X : Extended;
begin
for j:=1 to nn do
begin
X := A[p,j];
A[p,j] := A[q,j];
A[q,j] := X;
end;
X := B[p];
B[p] := B[q];
B[q] := X;
end;
Procedure Divid(p : Byte; k : Extended);
var
j : Byte;
begin
for j := 1 to nn do
A[p, j] := A[p, j] / k;
B[p] := B[p] / k;
end;
Procedure Sub(p, q : Byte);
var
j : Byte;
begin
for j := 1 to nn do
A[p, j] := A[p, j] - A[q, j];
B[p] := B[p] - B[q];
end;
Function CheckNullStr(p : Byte): Boolean;
var
i : Byte;
k : Extended;
begin
k := 0;
for i := 1 to nn do
k:=k + Abs (A[p, i]);
CheckNullStr := ((k=0) and (B[p]=0));
end;
Function BadStr(p : Byte) : Boolean;
var
i : Byte;
k : Extended;
begin
k := 0;
for i := 1 to nn do
k:=k + Abs (A[p, i]);
BadStr := ((k=0) and (B[p]<>0));
end;
Procedure KillNullStr(r : Byte);
begin
Swap(r, m);
dec(m);
end;
Function IsZeroDiag : Boolean;
var
i : Byte;
q : Boolean;
begin
q := false;
for i := 1 to m + 1 do
q := q or (A[i,i] = 0) ;
IsZeroDiag := q;
end;
Procedure FindNZD ;
var
i : Byte;
begin
for i := k+1 to m + 1 do
if A[i,k] <> 0.0 then
begin
Swap(k, i);
end;
end;
Procedure SolveMattr;
var
i : Byte;
Label
Loop;
begin
k := 1;
Loop:
While k <= m+1 do {???}
begin
for i := 1 to m+1 do
if not BadStr(i) then
begin
if CheckNullStr(i) then
begin
KillNullStr(i);
goto loop;
end;
end
else {ViewResult(0)};
FindNZD;
p:=0;
for i := k to nn do
if A[k,i] = 0 then inc(p)
else Break;
for i := 1 to m + 1 do
if A[i,k+p]<>0 then
Divid(i, A[i,k+p]);
for i :=1 to m + 1 do
if (i<>(k+p)) and (A[i,k+p]<>0.0) then
Sub(i,k+p);
inc(k);
end;
for i:=1 to m + 1 do if A[i,i]<>0 then
Divid(i,A[i,i]);
i := 1;
While i <= m + 1 do
begin
{ if BadStr(i) then {ViewResult(0);}
if CheckNullStr(i) then KillNullStr(i);
inc(i);
end;
end;
Function Xdg(x, dg : Extended):Extended;
begin
Xdg := exp(dg*ln(x));
end;
Function Fi(x:Extended):Extended;
var
tmp : Extended;
k : Word;
begin
tmp := 0;
for k := 0 to m do tmp := tmp + B[k+1]*Xdg(x, 1.0*k);
Fi := tmp;
end;
Procedure CalledKeyProc;
begin
x := x0;
dx := (xn-x0)/1000.0;
MoveTo(20,Yb-Round((Fi(x0)-ymin)/(ymax-ymin)*450));
SetColor(LightRed);
While x <= xn do
begin
LineTo(20+Round((x-x0)/(xn-x0)*600),Yb-Round((Fi(x)-ymin)/(ymax-ymin)*450));
x := x + dx;
end;
end;
Function Max(Ar : array of Extended) : Extended;
var
Mx : Word;
j : Word;
begin
Mx := 0;
for j := 0 to ndot do
if Ar[Mx]Ar[j] then Mn := j;
Min := Ar[Mn];
end;
Procedure SetMattr;
var
i, k : Word;
tmp : Extended;
begin
for s := 0 to m do {по строкам}
begin
for k := 0 to m do {По элементам строки}
begin
tmp := 0.0;
for i := 0 to ndot do
tmp := tmp + Xdg(xxx[i], 1.0*(k+s));
A[s+1, k+1] := tmp;
end;
tmp := 0.0;
for i := 0 to ndot do
tmp := tmp + y[i]*Xdg(xxx[i], 1.0*s);
B[s+1] := tmp;
end;
end;
Procedure ScreenMattr;
var
i, j : Byte;
begin
Writeln;
for i := 1 to m+1 do
begin
{ for j := 1 to nn do
Write(A[i, j]:1 :3,' ');}
WriteLn('a',i-1,':| ',B[i]:1 :3);
end;
ReadKey;
end;
Function Dispersia : Extended;
var
i : Longint;
tmp : Extended;
begin
tmp := 0; {!!!!!!!!!!!!!}
for i := 0 to ndot do
tmp := tmp + sqr(y[i]-Fi(xxx[i]));
Dispersia := tmp/(ndot+1);
end;
begin
if ParamCount=0 then
begin
WriteLn('Введите имя файла, содержащего данные...');
WriteLn('Структура файла:');
WriteLn('{Число точек}');
WriteLn('{Степень полинома}');
WriteLn('... {координаты точек} ...');
Exit;
end;
Assign(ff, ParamStr(1));
Reset(ff);
ReadLn(ff, ndot);
if ndot > n then
begin
Close(ff);
WriteLn('Слишком много точек в файле!!!');
Halt;
end;
Randomize;
ClrScr;
ReadLn(ff, m);
ReadLn(ff, x0);
ReadLn(ff, xn);
nn := m+1;
dec(ndot);
for i := 0 to ndot do
begin
ReadLn(ff, xxx[i]);
ReadLn(ff, y[i]);
end;
Close(ff);
dx := (xn-x0)/ndot;
WriteLn('Писло точек ', ndot+1);
WriteLn('Степень полинома ', m);
ymax := Max(y);
ymin := Min(y);
x0 := Min(xxx);
xn := Max(xxx);
if (xn <= 0) or (x0 <= 0) then
begin
for i := 0 to ndot do
xxx[i] := xxx[i] + Abs(xn)+Abs(x0);
x0 := Min(xxx);
xn := Max(xxx);
end;
WriteLn('x0 ', x0);
WriteLn('xn ', xn);
if (x0 = 0) and (xn = 0) then
begin
WriteLn('Длина отрезка по оси абсцисс равна 0 !!!');
Halt;
end;
SetMattr;
SolveMattr;
ScreenMattr;
Dr := Detect; {Инициализируем графику}
InitGraph(Dr, R, '');
e := GraphResult;
if e <> grOk then {Если ошибка}
begin
WriteLn('Ошибка инициализации графики!');
WriteLn(GraphErrorMsg(e));
ReadKey;
Halt($FFFF);
end;
MoveTo(20+Round((xxx[0]-x0)/(xn-x0)*600), Yb-Round((y[0]-ymin)/(ymax-ymin)*450));
SetColor(LightGreen);
Circle(20+Round((xxx[0]-x0)/(xn-x0)*600), Yb-Round((y[0]-ymin)/(ymax-ymin)*450), 1);
for i := 1 to ndot do
begin
SetColor(Yellow);
LineTo(20+Round((xxx[i]-x0)/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450));
SetColor(LightGreen);
Circle(20+Round((xxx[i]-x0)/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450), 1);
end;
ReadKey;
CalledKeyProc;
ReadKey;
CloseGraph;
dx := (xn-x0)/ndot;
for i := 0 to ndot do
WriteLn(i, ' : ', x0+dx*i:1:3, ' -> ', sqr(Fi(x0+dx*i)-y[i]));
Write('Дисперсия равна: ');
WriteLn(Dispersia:1);
Write('Среднеквадратичное отклонение равно: ');
WriteLn(sqrt(Dispersia):1);
ReadKey;
end.program
