|
Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)
outtextxy(80,130,'Меню ввода параметров нахождения');
outtextxy(80,140,' интеграла');
outtextxy(80,180,' Ввести количество узлов(n)');
outtextxy(80,210,' Ввести приделы интегрирования');
outtextxy(80,240,' Ввести функцию');
outtextxy(80,270,' Считать интеграл');
outtextxy(80,300,' Выход ');
end
else
begin
outtextxy(80,130,'Menu of entering the parameters');
outtextxy(80,140,' of integral');
outtextxy(80,180,' Put in the number of units ');
outtextxy(80,210,' Enter the bounds of integral');
outtextxy(80,240,' Enter function');
outtextxy(80,270,' Count integral');
outtextxy(80,300,' Exit ');
end;
end;
until c=#13;
c:='t';
case (abs(x) mod 5) of
0:
begin
wwodn(ea,n);
end;
1:
wwodab(ea,a,b);
2:
begin
helpwin(ea);
setcolor(15);
setfillstyle(1,9);
bar(70,200,340,300);
rectangle(75,205,335,295);
rectangle(77,207,333,293);
if ea mod 2 =0 then
begin
outtextxy(86,227,'Введите функцию f(x):');
setcolor(14);
outtextxy(360,140,' В этом окне необходимо');
outtextxy(360,155,' ввести саму функцию.');
outtextxy(360,200,'Примечание: 1.данная программа ');
outtextxy(360,215,'распознает только ');
outtextxy(360,230,'элементарные функции.');
outtextxy(360,245,'(x,cos(x) и др.)');
outtextxy(360,260,’2.При неправильном вводе’);
outtextxy(360,275,’по умолчанию f(x)=x;’);
outtextxy(360,275,’3.Если после нажатия ENTER’);
outtextxy(360,275,’ничего не произошло, то
outtextxy(360,275,’занововведите функцию.’);
end
else
begin
outtextxy(86,227,'Enter function f(x):');
setcolor(14);
outtextxy(360,140,' In this window you have');
outtextxy(360,155,' to enter the function.');
outtextxy(360,200,'Note: This version of ');
outtextxy(360,215,'programm can indentify only ');
outtextxy(360,230,'simple functions, as');
outtextxy(360,245,'x,cos(x) and other.');
end;
setfillstyle(1,0);
bar(86,255,330,275);
readln;
gotoxy(13,17);
read(st);
writeln(st);
readln;
end;
3:if (n0)and(a<>b)and(st<>'')and((abs(x) mod 5)=3);
end;
procedure win3(ea:word;n:integer;a,b:real;int:double;f:string;h:array of
double;var k:word);
{Последнее окно просмотра результатов}
var
i:integer;
c:char;
x:longint;
p1,p:string;
y:array[0..16] of double;
begin
funktia(n,a,b,y,1,f);
f:='('+f+')'+'dx =';
repeat
x:=-600000;
newsc(ea);
setfillstyle(1,2);
bar(170,120,490,360);
setcolor(14);
rectangle(175,125,485,355);
rectangle(177,127,483,353);
settextstyle(0,0,0);
setfillstyle(1,1);
bar(180,170,480,190);
if ea mod 2 =0 then
begin
outtextxy(180,135,Функция распознана.Интеграл подсчитан.');
outtextxy(180,180,' Посмотреть значение интеграла');
outtextxy(180,210,'Посмотреть коэффициенты Ньютона-Котеса');
outtextxy(180,240,' Посмотреть значения функции');
outtextxy(180,270,' Посмотреть график' );
outtextxy(180,300,' Считать снова');
outtextxy(180,330,' Выход ');
end
else
begin
outtextxy(180,135,'Function Indentified.Integral counted.');
outtextxy(180,180,' View value of integral');
outtextxy(180,210,' View Newton-Cotes coefficients');
outtextxy(180,240,' Veiw values of function');
outtextxy(180,270,' View graphik ' );
outtextxy(180,300,' Count again');
outtextxy(180,330,' Exit ');
end;
repeat
if keypressed then
begin
c:=readkey;
case c of
#80:
x:=x-1;
#72:
x:=x+1;
end;
setfillstyle(1,1);
case (abs(x) mod 6) of
0:
begin
bar(180,170,480,190);
setfillstyle(1,2);
bar(180,200,480,220);
bar(180,320,480,340);
end;
1:
begin
bar(180,200,480,220);
setfillstyle(1,2);
bar(180,170,480,190);
bar(180,230,480,250);
end;
2:
begin
bar(180,230,480,250);
setfillstyle(1,2);
bar(180,200,480,220);
bar(180,260,480,280);
end;
3:
begin
bar(180,260,480,280);
setfillstyle(1,2);
bar(180,230,480,250);
bar(180,290,480,310);
end;
4:
begin
bar(180,290,480,310);
setfillstyle(1,2);
bar(180,260,480,280);
bar(180,320,480,340);
end;
5:
begin
bar(180,320,480,340);
setfillstyle(1,2);
bar(180,290,480,310);
bar(180,170,480,190);
end;
end;
if ea mod 2 =0 then
begin
outtextxy(180,135,'Функция распознана.Интеграл подсчитан.');
outtextxy(180,180,' Посмотреть значение интеграла');
outtextxy(180,210,'Посмотреть коэффициенты Ньютона-Котеса');
outtextxy(180,240,' Посмотреть значения функции');
outtextxy(180,270,' Посмотреть график ' );
outtextxy(180,300,' Считать снова');
outtextxy(180,330,' Выход ');
end
else
begin
outtextxy(180,135,'Function Indentified.Integral counted.');
outtextxy(180,180,' View value of integral');
outtextxy(180,210,' View Newton-Cotes coefficients');
outtextxy(180,240,' Veiw values of function');
outtextxy(180,270,' View graphik ' );
outtextxy(180,300,' Count again');
outtextxy(180,330,' Exit ');
end;
end;
until c=#13;
c:='t';
case (abs(x) mod 6) of
0:begin
setcolor(15);
setfillstyle(1,12);
bar(140,200,490,280);
rectangle(145,205,485,275);
rectangle(147,207,483,273);
settextstyle(2,0,1);
setusercharsize(1,1,5,1);
outtextxy(170,210,'S');
settextstyle(2,0,4);
str(a:3:3,p);
outtextxy(160,257,p);
str(b:3:3,p);
outtextxy(160,212,p);
settextstyle(3,0,2);
outtextxy(180,224,f);
p:='';
str(abs(int):7:3,p);
outtextxy(190+length(f)*12,224,p);
readln;
end;
1:
begin
newsc(ea);
setfillstyle(1,2);
bar(170,120,490,180+n*15);
setcolor(14);
rectangle(175,125,485,175+n*15);
rectangle(177,127,483,173+n*15);
if ea mod 2 =0 then
begin
outtextxy(180,130,'Коэффициенты Ньютона-Котеса:');
outtextxy(180,140+(n+1)*15,'Нажмите ENTER для продолжения');
end
else
begin
outtextxy(180,130,'Newton-Cotes coefficients:');
outtextxy(180,140+(n+1)*15,'Press ENTER to continue');
end;
hkoef(n,h);
for i:=0 to n do
begin
str(i,p);str(h[i]:2:4,p1);
p:='H'+p+' = '+p1;
outtextxy(180,140+i*15,p);
end;
readln;
end;
2:begin
newsc(ea);
setfillstyle(1,2);
bar(170,120,490,180+n*15);
setcolor(14);
rectangle(175,125,485,175+n*15);
rectangle(177,127,483,173+n*15);
if ea mod 2 =0 then
begin
outtextxy(180,130,'Значения функции:');
outtextxy(180,140+(n+1)*15,'Нажмите ENTER для продолжения');
end
else
begin
outtextxy(180,130,'Values of function:');
outtextxy(180,140+(n+1)*15,'Press ENTER to continue');
end;
for i:=0 to n do
begin
str(i,p);str(y[i]:2:4,p1);
p:='Y'+p+' = '+p1;
p1:='';
outtextxy(180,140+i*15,p);
str((a+i*(b-a)/n):2:4,p1);
str(i,p);
if ea mod 2 = 0 then
p:=',При '+'X'+p+' = '+p1
else
p:=',When '+'X'+p+' = '+p1;
outtextxy(285,140+i*15,p);
end;
readln;
end;
3:
graphik(ea,a,b,f);
5:
begin
closegraph;
halt;
end;
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)71) then
line(320+i*24,240,320+i*24,242);
if ((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 ((k0)or(k10)) 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*c560 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)0) then
begin
if (abs(240-k1)>2) then
begin
if k17 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)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)71) then
line(320+i*24,20,320+i*24,460);
if ((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-ков на коэффициенты}
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;
end;
until c=#13;
repeat
newsc(ea);
win2(ea,k); {Ввод способа задания функции}
case k of
0:
wwod1(ea,y,n,a,b);
1:
begin
wwod2(ea,ea,n1,a1,b1,f);
n:=n1;a:=a1;b:=b1;
k:=4;
end;
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.
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
Рассмотрим результаты тестовых испытаний для функций 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) наблюдался обратный эффект.
Рекомендуемый диапозон n: от 7 до 13.
1) Интерфейс программы составлен на 2 языках: русском и английском. Переход
с одного языка на другой осуществляется в вводном окне путем нажатия
клавишы F2. Сменить язык можно только в этой части программы.
2) При вводе значений функции вручную необходимо вводить только цифры и
после каждого ввода нажимать клавишу ENTER.
3) При испытании программы под разные операционные системы(Dos, Windows 98-
2k,NT, из под паскаля) происходил непонятный баг с неверным выводом на
экран значений коэффициентов Ньютона-Котеса, хотя интеграл считался
верно. Для нормального нахождения их желательно запускать программу через
Dos.
4) При вводе параметров в “Меню задания параметров нахождения интеграла”
желательно их вводить постепенно сверху вниз, т.е. сначала ввести
количество узлов интерполяции, затем пределы интегрирования, а уж потом
вводить саму функцию.
5) Данная версия программы не способна распознавать все функции. Она может
распознать только стандартные функции Турбо Паскаля и еще несколько
дотполнительных: sin(x)/x, cos(x)*x ,arcctg(x). Для работы со
специфическими функциями необходимо в модуле K-unit в процедуре RASPOSN в
конце, перед end else, добавить :
k:=pos(‘Формула f(x)’,f);
if k<>0 then t:= ‘Формула f(x)’;
где ‘Формула f(x)’ – желаемая формула для распознования.
6) Вся помощь по вводу и работе с пограммой выводится в окне помощи.
Для нахождения интеграла существует много методов, однако, метод
Ньютона-Котеса один из самых быстрых: достаточно знать значения
коэффициентов для n=4, чтобы с точностью до сотых мгновенно посчитать
интеграл. Быстрота и простота –главные части этого метода.
В.И. Грызлов «Турбо Паскаль 7.0» Москва: ДМК 2000г.
Данилина «Численные методы» Москва: Высшая школа 1978г.
-----------------------
[pic]
[pic]
Страницы: 1, 2, 3
|
|