Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)
Рефераты >> Кибернетика >> Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)

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-ков на коэффициенты}


Страница: