Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)Рефераты >> Кибернетика >> Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)
funktia(1,x-1,x,y,c,f1);
k:=round(240-(y[0])*c);
k1:=round(240-(y[1])*c);
if ((k<480)and(k>0)or(k1<480)and(k1>0)) then
line(319-round((240+x1)/10)+x,k,320-round((240+x1)/10)+x,k1);
end;
if (v mod 2)=0 then
begin
funktia(1,a,b,y,1,f1);
k:=round(240-(y[0])*c);
k1:=round(240-(y[1])*c);
line(320-round((240+x1)/10)+round(a*c),k,320-round((240+x1)/10)+round(a*c),240);
line(320-round((240+x1)/10)+round(b*c),k1,320-round((240+x1)/10)+round(b*c),240);
if 320-round((240+x1)/10)+a*c<80 then
begin
funktia(1,-240/c,240/c,y,1,f1);
k:=round(240-(y[0])*c);
line(80,k,80,240);
end;
if 320-round((240+x1)/10)+b*c>560 then
begin
funktia(1,(-240-round((240+x1)/10))/c,(240-round((240+x1)/10))/c,y,1,f1);
k1:=round(240-(y[1])*c);
line(560,k1,560,240);
end;
for x:= -240 to 240 do
begin
funktia(1,x-1,x,y,c,f1);
k1:=round(240-(y[1])*c);
if ((x/c)>a) and ((x/c)<b) and(k1<>0) then
begin
if (abs(240-k1)>2) then
begin
if k1<240 then
k1:=k1+1
else
k1:=k1-1;
if c>7 then
setfillstyle(6,3)
else
setfillstyle(1,3);
floodfill(320-round((240+x1)/10)+x,k1,15);
end;
end;
end;
end;
str(x1,f2);
outtextxy(1,450,f2);
if (n mod 2)=0 then
for i:=-9 to 10 do
begin
settextstyle(2,0,2);
setcolor(14);
if ((320+i*24)<561) and ((320+i*24)>71)and(i<>0) then
begin
str((i*24+round((240+x1)/10))/c:2:2,f);
p:=247;
outtextxy(310+i*24,p,f);
str(-i*24/c:2:2,f);
outtextxy(330,240+i*24,f);
end;
end;
for i:=-9 to 10 do
begin
setcolor(15);
if ((r mod 2)=1) and (i<>0) then
begin
if ((320+i*24)<561) and ((320+i*24)>71) then
line(320+i*24,20,320+i*24,460);
if ((240+i*24)<461)and(240+i*24>19) then
line(80,240+i*24,560,240+i*24);
end;
end;
setcolor(15);
d:=readkey;
case d of
#75:
begin
x1:=x1-30;
x2:=x2-30;
end;
#77:
begin
x1:=x1+30;
x2:=x2+30;
end;
#80:
if c>1 then
c:=c-1;
#72:
c:=c+1;
#71:
c:=24;
#79:
n:=n+1;
#83:
r:=r+1;
#82:
v:=v+1;
#73:
begin
c:=24;
n:=0;r:=0;v:=0;x1:=-240;x2:=240;
end;
#59:
hwg(ea);
end;
until d=#13;
end;
end.
================================================
==========МОДУЛЬ UNIT==========
================================================
{$N+}
Unit k_unit;
{Модуль нахождения интеграл от многочлена q(q-1) (q-i+1)(q-i-1) (q-n),}
{где n-точность интеграла ,i-номер коофициента. }
interface
procedure rasposn(f:string;x:real;var ec:word;var t:real);
procedure hkoef(n:integer;var h:array of double);
procedure funktia(n:integer;a,b:real;var y:array of double;c:real;f:string);
procedure koef(w:array of double;n:integer;var e:array of double);
procedure mnogochlen(n,i:integer;var c:array of double);
function facktorial(n:integer):double;
function integral(w:array of double;n:integer):double;
function mainint(n:integer;a,b:real;y:array of double):double;
implementation
procedure rasposn(f:string;x:real;var ec:word;var t:real);
{Процедура распознования функции}
var
k:word;
begin
k:=pos('x',f);
if k<>0 then
begin {Распознавание функции}
ec:=1; {Код ошибки}
t:=x;
k:=pos('abs(x)',f);
if k<>0 then t:=abs(x);
k:=pos('sin(x)',f);
if k<>0 then t:=sin(x);
k:=pos('cos(x)',f);
if k<>0 then t:=cos(x);
k:=pos('arctg(x)',f);
if k<>0 then t:=arctan(x);
k:=pos('sqr(x)',f);
if k<>0 then t:=x*x;
k:=pos('exp(x)',f);
if k<>0 then t:=exp(x);
k:=pos('cos(x)*x',f);
if k<>0 then t:=cos(x)*x;
k:=pos('ln(x)',f);
if k<>0 then
begin
if x>0 then t:=ln(x)
else
t:=0;
end;
k:=pos('sqrt(x)',f);
if k<>0 then
if x>=0 then t:=sqrt(x)
else t:=0;
k:=pos('arcctg(x)',f);
if k<>0 then t:=pi/2-arctan(x);
k:=pos('sin(x)/x',f);
if k<>0 then if x<>0 then t:=sin(x)/x;
end
else
ec:=0;
end;
procedure funktia(n:integer;a,b:real;var y:array of double;c:real;f:string);
{Процедур подсчет Y-ков и распознавания функции}
var
t,h,x:real;
k,i:integer;
es:word;
begin
h:=(b-a)/n;
for i:=0 to n do
begin
x:=(a+h*i)/c;
rasposn(f,x,es,t);
y[i]:=t;
end;
end;
procedure koef(w:array of double;n:integer;var e:array of double);
{Изменение коофициентов для интеграла}
var
t:integer;
begin
for t:=1 to n do
e[t]:=w[t]/(n-t+2);
end;
procedure mnogochlen(n,i:integer;var c:array of double);
{процедура нахождения коофициентов при Q^n(q в степени n )}
var
k,j:integer;
d:array[1 100] of double;
begin
d[1]:=1;
for j:=1 to n do
begin {Вычисление коэффициентов при раскрытии q*(q-1)*(q-2)* *(q-n)}
d[j+1]:=d[j]*j*(-1);
if j>1 then
for k:=j downto 2 do
d[k]:=d[k]+d[k-1]*j*(-1);
end;
c[1]:=d[1]; {Деление многочлена на (q-i) по схеме Горнера}
for j:=1 to n+1 do
c[j]:=i*c[j-1]+d[j];
koef(c,n,c); {Изменение коэффициентов при интегрировании}
end;
function facktorial(n:integer):double;
{функция нахождения факториала }
var
t:integer;
s:double;
begin
s:=1;
if n=0 then
s:=1
else
for t:=1 to n do
s:=s*t;
facktorial:=s;
end;
function integral(w:array of double;n:integer):double;
{функция подсчета самого интеграла}
var
t,p:integer;
s,c:double;
begin
s:=0;p:=n;
for t:=0 to p+1 do
s:=s+w[t]*exp((p-t+2)*ln(p)); {Подсчет интеграла}
integral:=s;
end;
procedure hkoef(n:integer;var h:array of double);
{Процедура подсчета коэф. Ньютона-Котеса}
var
p,j,d,c,i:integer;
kq:array[0 20] of double;
s:array[0 20] of double;
begin
p:=n;
if (p mod 2)=1 then {Вычисление половины от всех вычислений коэффициентов}
d:=round((p-1)*0.5)
else
d:=round(0.5*p);
for i:=0 to n do
begin
mnogochlen(p,i,kq);
s[i]:=integral(kq,p); {Формирование массива из интегралов}
end;
for i:=0 to d do
begin
if ((p-i) mod 2) = 0 then
c:=1
else
c:=(-1);
h[i]:=(c*s[i])/(facktorial(i)*facktorial(p-i)*p);
h[p-i]:=h[i];
end;
end;
function mainint(n:integer;a,b:real;y:array of double):double;
{функция подсчета основного интеграла}
var
sum:double;
p,i:integer;
kq,h:array[0 20] of double;
begin
p:=n;
hkoef(n,h);
sum:=0;
for i:=0 to p do
sum:=sum+h[i]*y[i]; {Сумма произведений y-ков на коэффициенты}