Программа аппроксимации табличной функций с помощью интерполяционного полинома Лагранжа
Usescrt,graphABC;
Constnmax = 100;
n=20;
Vark,m,i,j: integer;
ct,u:real;
x,y,r,pt:array[1..nmax] of real;
xn,xk,mx,my,f:real;
x0,y0,i1,rad,xc,yc:integer;
s,t:string;
bln: boolean;
begin
SetWindowSize(1200, 800);
y[1]:= -1.84; y[11]:=-1.33;
y[2]:=-1.98; y[12]:=-1.47;
y[3]:=-1.72; y[13]:=-1.50;
y[4]:=-1.58; y[14]:=-2.65;
y[5]:=-1.69; y[15]:=-1.65;
y[6]:=-1.59; y[16]:=-1.87;
y[7]:=-1.58; y[17]:=-1.61;
y[8]:=-1.64; y[18]:=-1.86;
y[9]:=-1.55; y[19]:=-1.84;
y[10]:=-1.35; y[20]:=-1.91;
for i:=1 to 20 do begin
x[i]:=i/10;
end;
writeln('массивузлов:');
for i:=1 to n do
writeln(x[i]:0:2,' ');
writeln;
readln;
writeln('массивзначений:');
for i:=1 to n do
writeln(y[i]:0:2,' ');
writeln;
readln;
for i:=1 to n do R[i]:=0;
for i:=1 to n do
begin
ct:=1;pt[1]:=1;k:=1;
for j:=1 to n do
ifi<>j then
begin
ct:=ct*(x[i]-x[j]);
pt[k+1]:=1;
if k>1 then
for m:=k downto 2 do
pt[m]:=pt[m-1]-pt[m]*x[j];
pt[1]:=-pt[1]*x[j];
k:=k+1;
end;
u:=y[i]/ct;
for m:=1 to n do
r[m]:=r[m]+pt[m]*u;
end;
writeln('массивкоэффициентов:');
for i:=n downto 1 do
begin
writeln(r[i]:0:2,' ');
end;
readln;
clrscr;
x0:=30;{центрэкрана}
y0:=windowheight div 2;
xk:=3;{интервалпоХ}
mx:=(x0+1000)/xk;{масштаб по Х}
my:=y0/7;{по У}
line(0,y0,windowwidth,y0);{оси}
line(x0,0,X0,windowheight);
fori1:=1 to 10 do{максимальное количество засечек в одну сторону}
begin
line(x0+round(i1*mx),y0-3,x0+round(i1*mx),Y0+3); {засечкинаоси Х}
line(x0-round(i1*mx),y0-3,x0-round(i1*mx),Y0+3);
line(x0+3,y0-round(i1*my),x0-3,y0-round(i1*my)); {засечкинаоси Y}
line(x0+3,y0+round(i1*my),x0-3,y0+round(i1*my));
str(i1,s);
{подписьоси Х}
textout(x0+round(i1*mx),y0+10,s);
textout(x0-round(i1*mx),y0+10,'-'+s);
{подписьоси Y}
textout(x0-25,y0-round(i1*my),s);
textout(x0-25,y0+round(i1*my),'-'+s);
end;
{центр}
textout(x0+5,y0+10,'0');
{подписиконцовосей}
textout(windowwidth-10,y0-10,'X');
textout(x0+5,10, 'Y');
{график}
for i:=1 to n do
begin
str(i,t);
rad:= 5;
for xc:= x0+round(x[i]*mx) - rad to x0+round(x[i]*mx) + rad do
foryc:= y0-round(y[i]*my) - rad to y0-round(y[i]*my) + rad do begin
bln:= sqr(xc - (x0+round(x[i]*mx))) + sqr(yc - (y0-round(y[i]*my))) <= sqr(rad);
ifbln then SetPixel(xc, yc, clRed);
end;
TextOut(x0+round(x[i]*mx)+7,y0-round(y[i]*my)-5,t);
end;
end.
Результат выполнения программы:
Программа аппроксимации табличной функции кубическими сплайнами
Программа кубических сплайнов
usescrt,graphABC;
const n=20;
x0=0.1;
x9=2;
h=0.05;
typevec = array [0..100] of real;
vari,k:integer; x1,p,p1,p2:real;
x,f,c:vec;
procedure tab(n:integer; varx,f:vec);
var i:integer;
begin
for i:=1 to n do begin
x[i]:=i/10;
end;
SetWindowSize(1200, 800);
f[1]:= -1.84; f[11]:=-1.33;
f[2]:=-1.98; f[12]:=-1.47;
f[3]:=-1.72; f[13]:=-1.50;
f[4]:=-1.58; f[14]:=-2.65;
f[5]:=-1.69; f[15]:=-1.65;
f[6]:=-1.59; f[16]:=-1.87;
f[7]:=-1.58; f[17]:=-1.61;
f[8]:=-1.64; f[18]:=-1.86;
f[9]:=-1.55; f[19]:=-1.84;
f[10]:=-1.35; f[20]:=-1.91;
end;
procedurecs(n:integer; varx,f,c:vec);
vari,j,m:integer; a,b,r:real; k:vec;
begin k[1]:=0; c[1]:=0;
for i:=2 to n do begin j:=i-1; m:=j-1;
a:=x[i]-x[j]; b:=x[j]-x[m]; r:=2*(a+b)-b*c[j]; c[i]:=a/r;
k[i]:=(3*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r
end;
c[n]:=k[n];
for i:=n-1 downto 2 do c[i]:=k[i]-c[i]*c[i+1]
end;
procedure sp(n:integer; varx,f,c:vec; var x1,p,p1,p2:real);
vari,j:integer; a,b,d,q,r:real;
begin
i:=1; while (x1>x[i]) and (i<n) do i:=i+1;
j:=i-1; a:=f[j]; b:=x[j]; q:=x[i]-b;
r:=x1-b; p:=c[i]; d:=c[i+1];
b:=(f[i]-a)/q-(d+2*p)*q/3; d:=(d-p)/q*r;
p1:=b+r*(2*p+d); p2:=2*(p+d); p:=a+r*(b+r*(p+d/3))
end;
begin tab(n,x,f);
cs(n,x,f,c); k:=round((x9-x0)/h+1); x1:=x0;
writeln(' x ',' p ');
for i:=1 to k do begin sp(n,x,f,c,x1,p,p1,p2);
writeln(x1,' ',p:2:2); x1:=x1+h;
readln;
end;
end.
Результат выполнения программы:
Программа построения графика функции:
usescrt,graphABC;
const n=20;
x0=0.1;
x9=2;
h=0.05;//изменяя шаг, изменится число точек разбиения
typevec = array [0..100] of real;
vari,k:integer; x1,p,p1,p2:real;
xn,xk,mx,my,ct,u:real;
x00,y00,i1,rad,xc,yc:integer;
s,t:string;
bln: boolean;
x,f,c:vec;
procedure tab(n:integer; varx,f:vec);
var i:integer;
begin
for i:=1 to n do begin
x[i]:=i/10;
end;
SetWindowSize(1200, 800);
f[1]:= -1.84; f[11]:=-1.33;
f[2]:=-1.98; f[12]:=-1.47;
f[3]:=-1.72; f[13]:=-1.50;
f[4]:=-1.58; f[14]:=-2.65;
f[5]:=-1.69; f[15]:=-1.65;
f[6]:=-1.59; f[16]:=-1.87;
f[7]:=-1.58; f[17]:=-1.61;
f[8]:=-1.64; f[18]:=-1.86;
f[9]:=-1.55; f[19]:=-1.84;
f[10]:=-1.35; f[20]:=-1.91;
clrscr;
x00:=30;{центрэкрана}
y00:=windowheight div 2;
xk:=3;{интервал по Х}
mx:=(x00+1000)/xk;{масштаб по Х}
my:=y00/7;{поУ}
line(0,y00,windowwidth,y00);{оси}
line(x00,0,X00,windowheight);
for i1:=1 to 10 do{максимальное количество засечек в одну сторону}
begin
line(x00+round(i1*mx),y00-3,x00+round(i1*mx),Y00+3); {засечкинаосиХ}
line(x00-round(i1*mx),y00-3,x00-round(i1*mx),Y00+3);
line(x00+3,y00-round(i1*my),x00-3,y00-round(i1*my)); {засечкинаоси Y}
line(x00+3,y00+round(i1*my),x00-3,y00+round(i1*my));
str(i1,s);
{подписьосиХ}
textout(x00+round(i1*mx),y00+10,s);
textout(x00-round(i1*mx),y00+10,'-'+s);
{подписьоси Y}
textout(x00-25,y00-round(i1*my),s);
textout(x00-25,y00+round(i1*my),'-'+s);
end;
{центр}
textout(x00+5,y00+10,'0');
{подписиконцовосей}
textout(windowwidth-10,y00-10,'X');
textout(x00+5,10, 'Y');
end;
procedurecs(n:integer; varx,f,c:vec);
vari,j,m:integer; a,b,r:real; k:vec;
begin k[1]:=0; c[1]:=0;
for i:=2 to n do begin j:=i-1; m:=j-1;
a:=x[i]-x[j]; b:=x[j]-x[m]; r:=2*(a+b)-b*c[j]; c[i]:=a/r;
k[i]:=(3*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r
end;
c[n]:=k[n];
for i:=n-1 downto 2 do c[i]:=k[i]-c[i]*c[i+1]
end;
procedure sp(n:integer; varx,f,c:vec; var x1,p,p1,p2:real);
vari,j:integer; a,b,d,q,r:real;
begin
i:=1; while (x1>x[i]) and (i<n) do i:=i+1;
j:=i-1; a:=f[j]; b:=x[j]; q:=x[i]-b;
r:=x1-b; p:=c[i]; d:=c[i+1];
b:=(f[i]-a)/q-(d+2*p)*q/3; d:=(d-p)/q*r;
p1:=b+r*(2*p+d); p2:=2*(p+d); p:=a+r*(b+r*(p+d/3))
end;
begin tab(n,x,f);
cs(n,x,f,c); k:=round((x9-x0)/h+1); x1:=x0;
for i:=1 to k do begin sp(n,x,f,c,x1,p,p1,p2); x1:=x1+h;
begin
{график}
str(i,t);
rad:= 5;
for xc:= x00+round(x1*mx) - rad to x00+round(x1*mx) + rad do
foryc:= y00-round(p*my) - rad to y00-round(p*my) + rad do begin
bln:= sqr(xc - (x00+round(x1*mx))) + sqr(yc - (y00-round(p*my))) <= sqr(rad);
ifbln then SetPixel(xc, yc, clRed);
end;
TextOut(x00+round(x1*mx)+7,y00-round(p*my)-5,t);
end;
end;
end.
Результаты работы программы:
