Рефераты. Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале) p> 1: begin bar(80,200,330,220); setfillstyle(1,15); bar(80,170,330,190); bar(80,230,330,250); helpwin(ea); if ea mod 2 =0 then begin outtextxy(360,140,' Нажмите ENTER для ввода'); outtextxy(360,155,'приделов интегрирования.'); end else begin outtextxy(360,140,' Press ENTER to put in'); outtextxy(360,155,'the bounds of integral.'); end; end;

2: begin bar(80,230,330,250); setfillstyle(1,15); bar(80,200,330,220); bar(80,260,330,280); helpwin(ea); if ea mod 2 =0 then begin outtextxy(360,140,' Нажмите ENTER для ввода'); outtextxy(360,155,'функции.'); end else begin outtextxy(360,140,' Press ENTER to enter'); outtextxy(360,155,'function.'); end; end;

3: begin bar(80,260,330,280); setfillstyle(1,15); bar(80,230,330,250); bar(80,290,330,310); helpwin(ea); if ea mod 2 =0 then begin outtextxy(360,140,' Нажмите ENTER для начала'); outtextxy(360,155,'подсчета самого интеграла.'); end else begin outtextxy(360,140,' Press ENTER to begin'); outtextxy(360,155,'integral calculations.'); end; end;

4: begin bar(80,290,330,310); setfillstyle(1,15); bar(80,260,330,280); bar(80,170,330,190); helpwin(ea); end; end; setcolor(12); if ea mod 2 =0 then begin 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(ab)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)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(i0) 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 (i0) 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 k0 then begin {Распознавание функции} ec:=1; {Код ошибки} t:=x; k:=pos('abs(x)',f); if k0 then t:=abs(x); k:=pos('sin(x)',f); if k0 then t:=sin(x); k:=pos('cos(x)',f); if k0 then t:=cos(x); k:=pos('arctg(x)',f); if k0 then t:=arctan(x); k:=pos('sqr(x)',f); if k0 then t:=x*x; k:=pos('exp(x)',f); if k0 then t:=exp(x); k:=pos('cos(x)*x',f); if k0 then t:=cos(x)*x; k:=pos('ln(x)',f); if k0 then begin if x>0 then t:=ln(x) else t:=0; end; k:=pos('sqrt(x)',f); if k0 then if x>=0 then t:=sqrt(x) else t:=0; k:=pos('arcctg(x)',f); if k0 then t:=pi/2-arctan(x); k:=pos('sin(x)/x',f); if k0 then if x0 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 k4; 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 k0 then t:= ‘Формула f(x)’; где ‘Формула f(x)’ – желаемая формула для распознования.
6) Вся помощь по вводу и работе с пограммой выводится в окне помощи.

Для нахождения интеграла существует много методов, однако, метод
Ньютона-Котеса один из самых быстрых: достаточно знать значения коэффициентов для n=4, чтобы с точностью до сотых мгновенно посчитать интеграл. Быстрота и простота –главные части этого метода.

В.И. Грызлов «Турбо Паскаль 7.0» Москва: ДМК 2000г.
Данилина «Численные методы» Москва: Высшая школа 1978г.

-----------------------
[pic]

[pic]



Страницы: 1, 2, 3



2012 © Все права защищены
При использовании материалов активная ссылка на источник обязательна.