Program _Approc_;{Аппроксимация функций методом наименьших квадратов}
{$E-,N+,G+,D+}
Uses Graph, CRT;
type
 Sto = array [ 1..70 ] of Extended;
const
  n = 10;
  x0 : Extended = 1.0;
  xn : Extended = 4.0;
  Yb : Integer = 464;
var
  m, nn, k, p, s : Byte;
  Dr, R, e : Integer;
  y : array [0..n] of Extended;
  dx, ymax, ymin, x : Extended;
  A          : array [ 1..70 ] of Sto;{Массив }
  B          : Sto ;
  i, j : Longint;
Function f(x, d, g, p, q, m, k, l : Extended):Extended;{Сама функция}
begin
  f := (d+exp((p*x+q)*ln(x)))*exp(m*x*x+k*x+l);
end;
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;
  While x <= xn do
  begin
    PutPixel(20+Round((x-x0)/(xn-x0)*600),Yb-Round((Fi(x)-ymin)/(ymax-ymin)*450), LightRed);
    x := x + dx;
  end;
end;
Function Max(Ar : array of Extended) : Extended;
var
  Mx : Word;
  j : Word;
begin
  Mx := 0;
  for j := 0 to n 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 n do
            tmp := tmp + Xdg(x0+dx*i, 1.0*(k+s));
          A[s+1, k+1] := tmp;
        end;
          tmp := 0.0;
          for i := 0 to n do
            tmp := tmp + y[i]*Xdg(x0+dx*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 n do
    tmp := tmp + sqr(y[i]-Fi(x0+dx*i));
  Dispersia := tmp/(n+1);
end;
begin
  Randomize;
  m := 15;
  nn := m+1;
  ClrScr;
  dx := (xn-x0)/n;
  for i := 0 to n do
{  y[i] := 50 + (2.0*Random - 1.0)+i*1.0;}
  y[i] := f(x0+dx*i, 0, 1.2, 0.002, -1.001, 0, -0.001, 0);
  ymax := Max(y);
  ymin := Min(y);
  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, Yb-Round((y[0]-ymin)/(ymax-ymin)*450));
  SetColor(LightGreen);
  Circle(20, Yb-Round((y[0]-ymin)/(ymax-ymin)*450), 1);
  for i := 1 to n do
  begin
    SetColor(Yellow);
    LineTo(20+Round(i*dx/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450));
    SetColor(LightGreen);
    Circle(20+Round(i*dx/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450), 1);
  end;
  ReadKey;
  CalledKeyProc;
  ReadKey;
  CloseGraph;
  dx := (xn-x0)/n;
  for i := 0 to n 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.

Назад