Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале) (стр. 4 из 5)

outtextxy(285,140+i*15,p);

end;

readln;

end;

graphik(ea,a,b,f);

begin

closegraph;

halt;

end;

until (abs(x) mod 6)=4;

k:=abs(x) mod 6;

end;

end.

================================================

========МОДУЛЬ GRAPHIC========

================================================

unit graphic;

interface

uses

k_unit,crt,graph;

procedure hwg(ea:word);

procedure graphik(ea:word;a,b:real;f1:string);

implementation

procedure hwg(ea:word);

{Процедура окна помощи при графике}

var

f:string;

begin

settextstyle(0,0,0);

setfillstyle(1,3);

bar(150,100,390,380);

setcolor(0);

rectangle(153,103,387,377);

rectangle(155,105,385,375);

setcolor(14);

if ea mod 2 =0 then

begin

outtextxy(160,115,' ОКНО ПОМОЩИ');

outtextxy(160,140,' Для работы с графиком');

outtextxy(160,155,' используйте клавиши:');

outtextxy(160,180,' PAGE UP-первоначальный');

outtextxy(160,195,' вид графика;');

outtextxy(160,210,' HOME-начальный масштаб;');

outtextxy(160,225,' INSERT-включить/выключеть');

outtextxy(160,240,' заливку области;');

outtextxy(160,255,' DELETE-включить/выключеть');

outtextxy(160,270,' сетку;');

outtextxy(160,285,' END-показать/убрать цифры');

outtextxy(160,300,' F1- Помощь;');

outtextxy(160,315,' Стрелки ВВЕРХ/ВНИЗ- ');

outtextxy(160,330,' увеличение/уменьшение');

outtextxy(160,345,' масштаб .');

outtextxy(160,360,'Для возрата нажмите ENTER.');

end

else

begin

outtextxy(160,115,' HELP WINDOW');

outtextxy(160,140,' For the work with graphic');

outtextxy(160,155,' use this keys:');

outtextxy(160,180,' PAGE UP-Primery form of');

outtextxy(160,195,' graphik;');

outtextxy(160,210,' HOME-Primery scale;');

outtextxy(160,225,' INSERT-Turn on/off inking');

outtextxy(160,240,' the field;');

outtextxy(160,255,' DELETE-Turn on/off the');

outtextxy(160,270,' net;');

outtextxy(160,285,' END-View/delete the figures');

outtextxy(160,300,' F1- Help;');

outtextxy(160,315,' Arrows UP/DOWN-Increase/ ');

outtextxy(160,330,' lower the scale;');

outtextxy(160,360,'Press ENTER to continue.');

end;

readln;

setcolor(15);

end;

procedure graphik(ea:word;a,b:real;f1:string);

{процедура построения графиков}

var

f,f2:string;

d:char;

i,v,r:integer;

x1,x2,n,p,x:integer;

c,k,k1:longint;

y:array[0..1] of double;

begin

x1:=-240;

x2:=240;

c:=24;

setcolor(15);

n:=0;v:=0;r:=0;

repeat

cleardevice;

settextstyle(0,0,0);

if ea mod 2 =0 then

begin

outtextxy(10,1,'Нажмите F1 для помощи');

str(c/24:2:2,f);

f:='Масштаб '+f+':1';

end

else

begin

outtextxy(10,1,'Press F1 for help');

str(c/24:2:2,f);

f:='Scale '+f+':1';

end;

outtextxy(200,1,f);

settextstyle(3,0,1);

outtextxy(307,10,'y');

outtextxy(574,235,'x');

outtextxy(310,240,'0');

setlinestyle(1,7,100);

line(70,240,580,240);

line(320,20,320,460);

line(320,20,315,25);

line(321,20,326,25);

line(580,239,575,244);

line(580,240,575,235);

line(70,239,580,239);

line(321,20,321,460);

for i:=-9 to 10 do

begin

if ((320+i*24)<561) and ((320+i*24)>71) then

line(320+i*24,240,320+i*24,242);

if ((240+i*24)<461)and(240+i*24>19) then

line(320,240+i*24,322,240+i*24);

end;

setcolor(15);

for x:= -240+round((240+x1)/10) to 240+round((240+x1)/10) do

begin

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;

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;

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;

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;

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;

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

mainint:=sum*(b-a);

end;

end.

================================================

=======ОСНОВНАЯ ПРОГРАММА=======

================================================

{$N+}

program Newton_Cotes_metod;{Программа нахождения определенного интеграла}

uses {методом Ньютона-Котеса }

k_unit,k_graph,graph,crt;

const

t=15;

var

c:char;

a1,b1,a,b:real;

n1,v,r,n:integer;

h,y:array[0..t] of double;

ea,k:word;

int:double;

f:string;

begin

ea:=10;

v:=detect;

initgraph(v,r,'');

cleardevice;

newsc(ea);

winwin1;

setcolor(15);

outtextxy(380,430,'Нажмите F2 для смены языка.');

repeat

win1(ea);

settextstyle(3,0,1);

outtextxy(178,340,'Press Enter...');

delay(13000);

bar(178,340,350,365);

delay(13000);

if keypressed then {Смена языка}

begin

c:=readkey;

if c=#60 then

begin

ea:=ea+1;

newsc(ea);

winwin1;

setcolor(15);

if ea mod 2 =0 then

outtextxy(380,430,'Нажмите F2 для смены языка.')

else

outtextxy(380,430,'Press F2 key to change language.');

end;

until c=#13;

repeat

newsc(ea);

win2(ea,k); {Ввод способа задания функции}

case k of

wwod1(ea,y,n,a,b);

begin

wwod2(ea,ea,n1,a1,b1,f);

n:=n1;a:=a1;b:=b1;

k:=4;

end;

if k=4 then

funktia(n,a,b,y,1,f);

int:=mainint(n,a,b,y); {Вычисление интеграла}

hkoef(n,h);

proline(ea);

win3(ea,n,a,b,int,f,h,k); {Последнее меню вывода результатов}

until k<>4;

closegraph;

end.

Рассмотрим результаты тестовых испытаний для функций sin(x) на интервале [-5;3] и exp(x) на интервале [2;8]

n=1	n=2	n=3	n=4	n=5	n=7
Sin(x)	4,040017	3,02112	0,087629	1,779012	1,537481	1,246
Exp(x)	8965,041	3581,999	3271,82	3002,908	2990,644	2974,322
N=9	n=12
1,273561	1,27366
2973,593	2973,569

Видно, что при увеличении числа узлов интерполяции точность растет, однако при больших n (n>15) наблюдался обратный эффект.