Program _Approc_;{Аппроксимация функций методом наименьших квадратов}
{$E-,N+,G+,D+}
Uses Graph, CRT;
type
Tkf = array [0..6] of Extended;
Real = Extended;
type
Sto = array [ 1..70 ] of Extended;
const
n = 10;
x0 : Extended = 1.0;
xn : Extended = 4.0;
Yb : Integer = 464;
Eps : Extended = 1.0E-14;
dtb : Extended = 1.0E-4;
HotKey : Boolean = False;
var
mm, kk : Byte;
Dr, R, e : Integer;
y : array [0..n] of Extended;
dx, ymax, ymin, x : Extended;
OldInt09h : Pointer;
Pkf : ^Tkf;
A : array [ 1..70 ] of Sto;{Массив }
B : Sto ;
var
i, j : Longint;
{ d, g, p, q, m, k, l : Extended;}
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 : Real;
begin
for j:=1 to n 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 : Real);
var
j : Byte;
begin
for j := 1 to n do
begin
A[p, j] := A[p, j] / k;
if Abs(A[p, j]) <1.0E-10 then {Иначе возможен неправильный результат}
A[p, j] := 0.0;
end;
B[p] := B[p] / k;
if Abs(B[p]) <1.0E-10 then
B[p] := 0.0;
end;
Procedure Sub(p, q : Byte);
var
j : Byte;
begin
for j := 1 to n do
begin
A[p, j] := A[p, j] - A[q, j];
if Abs(A[p, j]) <1.0E-10 then {Иначе возможен неправильный результат}
A[p, j] := 0.0;
end;
B[p] := B[p] - B[q];
if Abs(B[p]) <1.0E-10 then
B[p] := 0.0;
end;
Function CheckNullStr(p : Byte): Boolean;
var
i : Byte;
k : Real;
begin
k := 0;
for i := 1 to n do
k:=k + Abs (A[p, i]);
CheckNullStr := ((k=0) and (B[p]=0));
end;
Function BadStr(p : Byte) : Boolean;
var
i : Byte;
k : Real;
begin
k := 0;
for i := 1 to n do
k:=k + Abs (A[p, i]);
BadStr := ((k=0) and (B[p]<>0));
end;
Procedure KillNullStr(r : Byte);
begin
Swap(r, mm);
dec(mm);
end;
Function IsZeroDiag : Boolean;
var
i : Byte;
q : Boolean;
begin
q := false;
for i := 1 to mm do
q := q or (A[i,i] = 0) ;
IsZeroDiag := q;
end;
Procedure FindNZD ;
var
i : Byte;
begin
for i := kk+1 to mm do
if A[i,kk] <> 0.0 then
begin
Swap(kk, i);
end;
end;
Procedure SolveMatr;
var
i : Byte;
Label
Loop{:Word};
begin
kk := 1;
Loop:
While kk <= mm do {???}
begin
for i := 1 to mm 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 := kk to n do
if A[k,i] = 0 then inc(p)
else Break;
for i := 1 to m do
if A[i,k+p]<>0 then
Divid(i, A[i,k+p]);
for i :=1 to m do
if (i<>(k+p)) and (A[i,k+p]<>0.0) then
Sub(i,k+p);
inc(kk);
end;
for i:=1 to m do if A[i,i]<>0 then
Divid(i,A[i,i]);
i :=1;
While i <= m do
begin
if BadStr(i) then ViewResult(0);
if CheckNullStr(i) then KillNullStr(i);
inc(i);
end;
end;
{Function S(d, g, p, q, m, k, l : Extended) : Extended; {Сумма квадратов отклонений}
{var
i : Longint;
Si : Extended;
begin
Si := 0.0;
for i := 0 to n do
Si := Si + sqr(f(x0+dx*i, d, g, p, q, m, k, l)-y[i]);
S := Si;
end;
{Procedure Approch(var d, g, p, q, m, k, l : Extended);
var
i : Longint;
t, dt : Extended;
d1, g1, p1, q1, m1, k1, l1 : Extended;
begin
dt := dtb;
t := l;
Write('l...');
if (S(d, g, p, q, m, k, t+dt)-S(d, g, p, q, m, k, t)>0.0) then dt := -dt;
While (S(d, g, p, q, m, k, t+dt)-S(d, g, p, q, m, k, t)<0.0) do
t := t + dt;
l1 := t;
WriteLn(t:5:5);
dt := dtb;
t := k;
Write('k...');
if (S(d, g, p, q, m, t+dt, l)-S(d, g, p, q, m, t, l)>0.0) then dt := -dt;
While (S(d, g, p, q, m, t+dt, l)-S(d, g, p, q, m, t, l)<0.0) do
t := t + dt;
k1 := t;
WriteLn(t:5:5);
dt := dtb;
t := m;
Write('m...');
if (S(d, g, p, q, t+dt, k, l)-S(d, g, p, q, t, k, l)>0.0) then dt := -dt;
While (S(d, g, p, q, t+dt, k, l)-S(d, g, p, q, t, k, l)<0.0) do
t := t + dt;
m1 := t;
WriteLn(t:5:5);
dt := dtb;
t := q;
Write('q...');
if (S(d, g, p, t+dt, m, k, l)-S(d, g, p, t, m, k, l)>0.0) then dt := -dt;
While (S(d, g, p, t+dt, m, k, l)-S(d, g, p, t, m, k, l)<0.0) do
t := t + dt;
q1 := t;
WriteLn(t:5:5);
dt := dtb;
t := p;
Write('p...');
if (S(d, g, t+dt, q, m, k, l)-S(d, g, t, q, m, k, l)>0.0) then dt := -dt;
While (S(d, g, t+dt, q, m, k, l)-S(d, g, t, q, m, k, l)<0.0) do
t := t + dt;
p1 := t;
WriteLn(t:5:5);
dt := dtb;
t := g;
WriteLn('g...');
if (S(d, t+dt, p, q, m, k, l)-S(d, t, p, q, m, k, l)>0.0) then dt := -dt;
While (S(d, t+dt, p, q, m, k, l)-S(d, t, p, q, m, k, l)<0.0) do
begin
t := t + dt;
GotoXY(1, WhereY);
Write(S(d, t+dt, p, q, m, k, l)-S(d, t, p, q, m, k, l) :5 :5, ' ', t :5 :5, ' ');
end;
g1 := t;
WriteLn(t:5:5);
dt := dtb;
t := d;
WriteLn('d...');
if (S(t+dt, g, p, q, m, k, l)-S(t, g, p, q, m, k, l)>0.0) then dt := -dt;
While (S(t+dt, g, p, q, m, k, l)-S(t, g, p, q, m, k, l)<0.0) do
begin
t := t + dt;
GotoXY(1, WhereY);
Write(S(t+dt, g, p, q, m, k, l)-S(t, g, p, q, m, k, l) :5 :5, ' ', t :5 :5, ' ');
end;
d1 := t;
WriteLn(t:5:5);
d := d1; g := g1; p := p1; q := q1; m := m1; k := k1; l := l1;
end;
Function Sgrad(k : array of Extended) : Extended; {Сумма квадратов отклонений}
{var
i : Word;
Si : Extended;
begin
Si := 0.0;
for i := 0 to n do
Si := Si + sqr(f(x0+dx*i, k[0], k[1], k[2], k[3], k[4], k[5], k[6])-y[i]);
Sgrad := Si;
end;
Function norma(k : array of Extended) : Extended;
var
i : Word;
Si : Extended;
begin
Si := 0;
for i := 0 to 6 do Si := sqr(k[i]);
norma := Si;
end;}
Procedure CalledKeyProc;FAR;
var
kkf : ^Tkf;
PPP : Pointer;
begin
PPP := Pkf;
kkf := PPP;
x := x0; dx := (xn-x0)/1000;
ymax := f(x0, kkf^[0],kkf^[1],kkf^[2],kkf^[3],kkf^[4],kkf^[5],kkf^[6]);
ymin := f(xn, kkf^[0],kkf^[1],kkf^[2],kkf^[3],kkf^[4],kkf^[5],kkf^[6]);
While x <= xn do
begin
PutPixel(20+Round((x-x0)/(xn-x0)*600),
Yb-Round((f(x,kkf^[0],kkf^[1],kkf^[2],kkf^[3],kkf^[4],kkf^[5],kkf^[6])-ymin)/(ymax-ymin)*450.0),
LightRed);
x := x + dx;
end;
{ SetColor(Yellow);
Circle(100, 100, 100);}
end;
{Procedure KeyProc;FAR;Assembler;
asm
pusha
push ds
push es
in al, 60h
cmp al, 1
jne @handler
mov byte ptr hotkey, al
@handler:
push offset @xfurther
push cs
jmp cs : dword ptr CalledKeyProc
@xfurther:
in al, 61h
push ax
or al, 80h
out 61h, al
pop ax
out 61h, al
mov al, 20h
out 20h, al
pop es
pop ds
popa
iret
end;}
{Procedure ApprochGrad(var d, g, p, q, m, k, l : Extended);
var
i : Word;
t, dt : Extended;
kf : array [0..6] of Extended;
gradf : array [0..6] of Extended;
Sg : Extended;
begin
Pkf := @kf;
kf[0]:=d;kf[1]:=g;kf[2]:=p;kf[3]:=q;
kf[4]:=m;kf[5]:=k;kf[6]:=l;
{ asm
push offset @further
push cs
jmp cs : dword ptr CalledKeyProc
@further:
end;
Exit;}
{ asm
push es
mov ax, 3509h
int 21h
mov word ptr oldint09h, bx
mov word ptr oldint09h+2, es
pop es
push ds
mov ax, seg KeyProc
mov ds, ax
mov dx, offset KeyProc
mov ax, 2509h
int 21h
pop ds
end;}
{ dt := 1.0E-8;
{ for i := 0 to 6 do kf[i] := 0;}
{ Repeat
for i := 0 to 6 do
begin
kf[i] := kf[i]+dt;
Sg := Sgrad(kf);
kf[i] := kf[i]-dt;
gradf[i] := (Sg-Sgrad(kf))/dt;
end;
Sg := norma(gradf);
for i := 0 to 6 do
begin
if gradf[i]>0.0 then kf[i] := kf[i]-dt*(1000/Sg+1);
if gradf[i]<0.0 then kf[i] := kf[i]+dt*(1000/Sg+1);
end;
{ GotoXY(1, WhereY);
Write(Sg:5:15);}
{ if HotKey then Break;}
{ if KeyPressed then if ReadKey=#27 then Break else CalledKeyProc;
Until Sg < Eps;
d := kf[0]; g := kf[1]; p := kf[2]; q := kf[3]; m := kf[4]; k := kf[5]; l := kf[6];
{ asm
push ds
mov ax, 2509h
mov dx, word ptr oldint09h+2
mov ds, dx
mov dx, word ptr oldint09h
pop ds
end;}
{ While KeyPressed do ReadKey;
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;
begin
{ WriteLn(Seg(ApprochGrad),':',Ofs(ApprochGrad));
Halt;}
{ d := 0; g := 1.2; p := 0.002; q := -1.001; m := 0; k := -0.001; l := 0;
d := 0; g := 0; p := 0; q := 0; m := 0; k := 0; l := 0;}
ClrScr;
dx := (xn-x0)/n;
for i := 0 to n do
y[i] := f(x0+dx*i, 0, 1.2, 0.002, -1.001, 0, -0.001, 0);
ymax := Max(y);
ymin := Min(y);
{ Approch(d, g, p, q, m, k, l);}
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), 2);
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), 2);
end;
ApprochGrad(d, g, p, q, m, k, l);
ReadKey;
CloseGraph;
{ WriteLn(d:5:5, ' > ', g:5:5, ' > ', p:5:5, ' > ', q:5:5, ' > ', m:5:5, ' > ', k:5:5, ' > ', l:5:5);}
end.
< -->