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




                    Министерство Высшего Образования РФ.


                   Московский Институт Электронной Техники
                          (Технический Университет)


                                 Лицей №1557



                               КУРСОВАЯ РАБОТА



                       “Вычисление интеграла  методом
                               Ньютона-Котеса”



                           Написал: Коноплев А.А.


                        Проверил: доцент Колдаев В.Д.



                               Москва, 2001г.



1.
   Введение..................................................................
   ................... 3
2.                                                             Теоретическая
   часть...................................................................4
3.                                                                  Алгоритм
   работы....................................................................
   ....8
4.                                                                       Код
   программы.................................................................
   ........17
 .                                                                    Модуль
   K_graph............................................................17
 .                                                                    Модуль
   Graphic.............................................................34
 .                                                                    Модуль
   K_unit...............................................................38
 . Основная программа....................................................40
5.                                                                  Тестовые
   испытания.................................................................
   42
6. Полезные советы по работе с программой.............................42
7.              Окна              ввода               и               вывода
   программы.............................................
8.
   Вывод.....................................................................
   .....................43
9.                                                                    Список
   литературы................................................................
   ...44



     Математика - одна из самых древних наук. Труды многих ученых вошли в
мировой фонд и стали основой современных алгебры и геометрии. В конце XVII
в., когда развитие науки шло быстрыми темпами, появились понятия
дифференцирование, а вслед за ним и интегрирование. Многие правила
нахождения неопределенного интеграла в то время не были известны, поэтому
ученые пытались  найти другие, обходные пути поиска значений. Первым
методом явился метод Ньютона – поиск интеграла через график функции, т.е.
нахождение площади под графиком, методом прямоугольников, в последствии
усовершенствованный в метод трапеций. Позже был придуман параболический
метод или метод Симпсона. Однако часть ученых терзал вопрос: А можно ли
объединить все эти методы в один??
Ответ на него был дан одновременно двумя математиками Ньютоном и Котесом.
Они вывели общую формулу, названную в их честь. Однако   их метод был
частично забыт. В этой работе будут изложены основные положения теории,
рассмотрены различные примеры,  приведены таблицы, полученные при различных
погрешностях, и конечно описана работа и код программы, рассчитывающей
интеграл методом Ньютона-Котеса.



Пусть некоторая функция f(x) задана в уздах интерполяции:
[pic]         (i=1,2,3…,n) на отрезке [а,b] таблицей значений:



|X0=a          |X1            |X2            |…             |XN=b          |
|Y0=f(x0)      |Y1=f(x1)      |Y2=f(x2)      |…             |YN=f(xN)      |

Требуется найти значение интеграла    [pic] .
  Для начала составим интерполяционный многочлен Лагранджа:

                                    [pic]

Для равноотстоящих узлов интерполяционный многочлен имеет вид:

                                    [pic]


где q=(x-x0)/h – шаг интерполяции, заменим подынтегральную функцию f(x)
интерполяционным многочленом Лагранжа:



                                    [pic]



Поменяем знак суммирования и интеграл и вынесем за знак интеграла
постоянные элементы:

                                    [pic]


      Так как dp=dx/h, то, заменив пределы интегрирования, имеем:


                                    [pic]

Для равноотстоящих узлов интерполяции на отрезке [a,b] величина  шаг
определяется как h=(a-b)/n. Представив это выражение для h в формулу (4) и
вынося (b-a) за знак суммы, получим:
                                    [pic]


Положим, что

                                    [pic]
где i=0,1,2…,n; Числа Hi называют коэффициентами Ньютона-Котеса. Эти
коэффиценты не зависят от вида f(x), а являются функцией только по n.
Поэтому их можно вычислить заранее. Окончательная формула выглядит так:
                                    [pic]
Теперь рассмотрим несколько примеров.



Пример 1.
   Вычислить с помощью метода Ньютона-Котаса:                  [pic]
  , при n=7.
        Вычисление.
   1) Определим шаг:  h=(7-0)/7=1.
   2)Найдем значения y:

|x0=0  |y0=1         |
|x1=1  |y1=0.5       |
|x2=2  |y2=0.2       |
|x3=3  |y3=0.1       |
|x4=4  |y4=0.0588    |
|x5=5  |y5=0.0384    |
|x6=6  |y6=0.0270    |
|x7=7  |y7=0.02      |



   3) Находим коэффициенты Ньютона-Котеса:
H1=H7=0.0435, H1=H6=0.2040, H2=H5=0.0760 ,H3=H4=0.1730
Подставим значения в формулу и получим:

При подсчете с помощью формулы Ньютона-Лейбница получим:



Пример 2.
Вычислить при помощи метода Ньютона-Котеса
[pic]     , взяв n=5;
Вычисление:
1) Определим  шаг h=(8-4)/5=0.8
2) Найдем значения y:

|x0=0     |y0=-2.61     |
|x1=4.8   |y1=0.42      |
|x2=5.6   |y2=4.34      |
|x3=6.4   |y3=6.35      |
|x4=7.2   |y4=4.38      |
|x5=8     |y5=-0.16     |


3) Находим коэффициенты Ньютона –Котеса:
H0=H5=0.065972 ;H1=H4=0.260417 ;H2=H3=0.173611 ;
4)Подставим значения в формулу и получим:
[pic]


    Рассмотрим частные случаи формулы Ньйтона-Котеса.
Пусть n=1 тогда
H0=H1=0.5 и конечная формула примет вид:
[pic] Тем самым в качестве частного случая нашей формулы мы получили
формулу трапеций.
Взяв n=3, мы получим
[pic] . Частный случай формулы Ньютона –Котеса – формула Симпсона



     Теперь произведем анализ алгоритма и рассмотрим  основной принцип
работы программы.
      Для вычисления интеграла сначала находятся коэффициенты Ньютона-
Котеса. Их нахождение осуществляется в процедуре hkoef.
Основной проблемой вычисления коэффициентов является интеграл от
произведения множителей. Для его расчета необходимо:

А) посчитать коэффициенты при раскрытии скобок при q
 (процедура mnogoclen)
Б) домножить их на 1/n , где n –степень при q (процедура koef)
В) подставить вместо q значение n (функция integral)

      Далее вычисляем факториалы (функция faktorial) и перемножаем
полученные выражения (функция mainint). Для увеличения быстроты работы
вводится вычисление половины от количества узлов интерполяции и последующей
подстановкой их вместо неподсчитанных.

               Процедура koef(w: массив;n:целый;var e:массив);

                                    [pic]

                   Процедура hkoef(n:целый;var h:массив);

                                    [pic]



                                    [pic]



                Процедура mnogochlen(n,i:целые;var c:массив );
                                    [pic]



                                    [pic]



   Процедура funktia(n:целая;a,b:вещест.;var y:массив;c:вещест.;f:строка);

                                    [pic]



                    Функция facktorial(n:целый):двойной;

                                    [pic]



                 Функция integral(w:массив;n:целый):двойной;

                                    [pic]



           Функция mainint(n:целый;a,b:вещест.;y:массив):двойной;

                                    [pic]



                             Основная программа

                                    [pic]



 Программа состоит из 8 файлов:
 . K_main.exe – файл загрузки основной программы
 . K_unit.tpu – модуль вычислительных процедур и функций
 . K_graph.tpu – модуль графических процедур
 . Graphic.tpu – модуль процедур для построения графика
 . Egavga.bgi – файл графической инициализации
 . Sans.chr, litt.chr – файлы шрифтов
 . Keyrus.com (не обязательно) – файл установки русского языка.
  Для работы программы с русским интерфайсом желательно запускать ее в
режиме DOS.

              ================================================
                      ==========МОДУЛЬ GRAPH==========
              ================================================
{$N+}
unit k_graph;
interface
uses
crt,graph,k_unit,graphic;
procedure winwin1;
procedure proline(ea:word);
procedure winwwodab(ea:word);
procedure error1(ea:word);
procedure helpwin(ea:word);
procedure error(ea:word);
procedure newsctext(ea:word);
procedure newsc(ea:word);
procedure win1(ea:word);
procedure win2(ea:word;var k:word);
procedure wwodn(ea:word;var n:integer);
procedure wwodab(ea:word;var a,b:real);
procedure wwod1(ea:word;var y:array of double;var n:integer;var a,b:real);
procedure wwod2(ea:word;var ea1:word;var n:integer;var a,b:real;var
st:string);
procedure win3(ea:word;n:integer;a,b:real;int:double;f:string;h:array of
double;var k:word);
implementation
procedure proline(ea:word);
{Проседура полосы процесса}
var
i:integer;
f:string;
c:char;
begin
     newsc(ea);
     setcolor(15);
     setfillstyle(1,7);
     bar(160,150,460,260);
     rectangle(165,155,455,255);
     rectangle(167,157,453,253);
     case (ea mod 2) of
          0: outtextxy(180,170,'     Идет работа .Ждите..');
          1: outtextxy(180,170,'    Working.Please wait..');
     end;
     setfillstyle(1,12);
     setcolor(0);
     rectangle(200,199,401,221);
     for i:=1 to 9 do
         line(200+i*20,200,200+i*20,220);
     delay(20000);
     for i:=1 to 100 do
     begin
         if ((i-1) mod 10)=0 then
            line(200+((i-1) div 10)*20,200,200+((i-1) div 10)*20,220);
         bar(round(200+2*(i-0.5)),200,200+2*i,220);
         delay(1100);
         setcolor(15);
         setfillstyle(1,7);
         bar(280,230,323,250);
         str(i,f);
         f:=f+'%';
         outtextxy(290,235,f);
         if (i mod 25) =0 then
            bar(170,180,452,198);
         if (ea mod 2)=0 then
         case (i div 25) of
              0:
                outtextxy(170,190,'Подготовка ');
              1:
                outtextxy(170,190,'Расчет коеффициентов в многочлене');
              2:
                outtextxy(170,190,'Расчет коеффициентов Ньютона-Котеса');
              3:
                outtextxy(170,190,'Расчет интеграла');
         end
         else
             case (i div 25) of
              0:
                outtextxy(170,190,'Prepearing');
              1:
                outtextxy(170,190,'Calculation of mnogochlen coeff.');
              2:
                outtextxy(170,190,'Calculation of Newton-Cotes coeff. ');
              3:
                outtextxy(170,190,'Calculation of integral');
                end;
         setfillstyle(1,12);
         setcolor(0);
     end;
end;
procedure winwwodn(ea:word);
{Окно ввода числа узлов интерполяции}
var
c:char;
f:string;
begin
     helpwin(ea);
     if (ea mod 2) =0 then
     begin
     outtextxy(360,140,'  В этом окне необходимо     ');
     outtextxy(360,155,' ввести количество узлов     ');
     outtextxy(360,170,' интерполяции, от которого   ');
     outtextxy(360,185,' будет зависить точность     ');
     outtextxy(360,200,' вычисления интеграл  и      ');
     outtextxy(360,215,' количество зн чений функции.');
     outtextxy(360,240,' ВНИМАНИЕ : НАСТОЯТЕЛЬНО     ');
     outtextxy(360,250,' РЕКОМЕНДУЕТСЯ НЕ ВВОДИТЬ    ');
     outtextxy(360,260,' ЗНАЧЕНИЕ N БОЛЬШЕ 12 !!     ');
     end
     else
     begin
     outtextxy(360,140,'  In this window you have to  ');
     outtextxy(360,155,' put into the number.         ');
     outtextxy(360,170,' The accuracy of calculation  ');
     outtextxy(360,185,' and the number of function   ');
     outtextxy(360,200,' parameters will depend on    ');
     outtextxy(360,215,' this number.                 ');
     outtextxy(360,240,' WARNING: IT IS HARDLY        ');
     outtextxy(360,250,' RECOMENDED NOT TO PUT IN     ');
     outtextxy(360,260,' NUMBER MORE THEN 12 !!       ');
     end;
     setcolor(2);
     setfillstyle(1,14);
     bar(70,200,340,300);
     rectangle(75,205,335,295);
     rectangle(77,207,333,293);
     if (ea mod 2) =0 then
     begin
     outtextxy(90,227,'Введите количество узлов(n):');
     outtextxy(80,270,'ВНИМАНИЕ: При больших n возможна');
     outtextxy(80,280,'некорректная работа компьютера!!');
     end
     else
     begin
     outtextxy(80,217,'Put in number of');
     outtextxy(80,227,'        interpolation units:');
     outtextxy(80,270,'WARNING:if you use big number  ');
     outtextxy(80,280,'of units,PC wont work properly!');
     end;
     setfillstyle(1,0);
     bar(190,240,230,255);
end;
procedure wwodn(ea:word;var n:integer);
{Процедура  ввода узлов n}
var
ec,p:integer;
k,f:string;
x:integer;
c:char;
begin
     newsc(ea);
     winwwodn(ea);
     repeat
     repeat
     winwwodn(ea);
     gotoxy(25,16);
     read(k);
     val(k,p,ec);
     if ec<>0 then
     begin
        error1(ea);
        readln;
     end;
     until ec=0;
     n:=p;
     if n>12 then
     begin
        if keypressed then
           c:=readkey;
        c:='r';
        setcolor(15);
        setfillstyle(1,12);
        bar(140,210,490,300);
        rectangle(145,215,485,295);
        rectangle(147,217,483,293);
        if (ea mod 2) =0 then
        begin
        outtextxy(150,227,'              Предупреждение!');
        outtextxy(150,237,'   Вы дейcтвительно хотите использовать');
        outtextxy(150,250,'          большое значение N ???');
        end
        else
            begin
                 outtextxy(150,227,'                Warning!!
');
                 outtextxy(150,237,'      Do you realy want to use a big
');
                 outtextxy(150,250,'     number interpolation units(N)???
');
                 end;
        sound(600);
        delay(4000);
        nosound;
        setfillstyle(1,2);
        bar(320,260,350,280);
        setfillstyle(1,12);
        bar(250,260,280,280);
        repeat
        if keypressed then
        begin
          c:=readkey;
          if (c=#80) or (c=#72) or (c=#77) or (c=#75) then
             x:=x+1;
          setfillstyle(1,2);
          if (x mod 2)=0 then
             begin
                  bar(250,260,280,280);
                  setfillstyle(1,12);
                  bar(320,260,350,280);
             end
          else
              begin
                   bar(320,260,350,280);
                   setfillstyle(1,12);
                   bar(250,260,280,280);
              END;

        end;
        if (ea mod 2) =0 then
        begin
        outtextxy(255,267,'ДА');
        outtextxy(325,267,'НЕТ');
        end
        else
        begin
             outtextxy(255,267,'YES');
             outtextxy(325,267,'NO');
        end;
        until c=#13;
        if abs(x mod 2)=1 then
        begin
               n:=0;
               setcolor(15);
               setfillstyle(1,2);
               bar(160,200,460,280);
               rectangle(165,205,455,275);
               rectangle(167,207,453,273);
               if (ea mod 2)=0 then
               begin
               outtextxy(180,227,'Для работы программы необходимо');
               outtextxy(180,237,'       заново ввести N.');
               outtextxy(180,247,' Нажмите ENTER для продолжения.');
               end
               else
               begin
                    outtextxy(180,227,' To continue you have to ');
                    outtextxy(180,237,'      again put in N.  ');
                    outtextxy(180,247,' Press ENTER to continue.');
               end;
               readln;
               readln;
        end;
     end;
     until n>0;
end;


procedure winwwodab(ea:word);
{Окно ввода приделов интегрирования}
var
f:string;
begin
     helpwin(ea);
     if (ea mod 2)=0 then
     begin
     outtextxy(360,140,'  В этом окне необходимо');
     outtextxy(360,155,' ввести сначала нижнее');
     outtextxy(360,170,' значение интеграл  и нажать');
     outtextxy(360,185,' ENTER, а затем ввести');
     outtextxy(360,200,' верхнее значение интеграла');
     outtextxy(360,215,' и снова нажать ENTER.');
     end
     else
     begin
     outtextxy(360,140,' In this window you have to:');
     outtextxy(360,155,'firstly, put in lower value  ');
     outtextxy(360,170,'of integral and press ENTER,');
     outtextxy(360,185,'then put in higher value');
     outtextxy(360,200,'of integral and press ENTER');
     end;
     setcolor(2);
     setfillstyle(1,5);
     bar(10,210,335,320);
     rectangle(15,215,330,315);
     rectangle(17,217,328,313);
     settextstyle(0,0,0);
     if (ea mod 2)=0 then
     begin
     outtextxy(20,230,'  Введите нижнее значение');
     outtextxy(20,244,' интеграл :');
     outtextxy(20,262,'  Введите верхнее значение');
     outtextxy(20,272,'интеграл :');
     end
     else
     begin
     outtextxy(20,230,'  Put in lower value of');
     outtextxy(20,244,' integral:');
     outtextxy(20,262,'  Put in higher value of');
     outtextxy(20,272,'integral:');
     end;
end;
procedure wwodab(ea:word;var a,b:real);
{Процедура  ввода  приделов интегрирования}
var
f:string;
k:string;
ec:integer;
begin
     newsc(ea);
     winwwodab(ea);
     readln;
     repeat
     winwwodab(ea);
     gotoxy(16,16);
     read(k);
     val(k,a,ec);
     if ec<>0 then
        error1(ea);
     until ec=0;
     readln;
     repeat
     winwwodab(ea);
     str(a:4:2,f);
     outtextxy(120,244,f);
     gotoxy(16,18);
     read(k);
     val(k,b,ec);
     if ec<>0 then
        error1(ea);
     until ec=0;
end;
procedure helpwin(ea:word);
{основа окна помощи}
begin
     setfillstyle(1,3);
     bar(350,100,590,380);
     setcolor(0);
     rectangle(353,103,587,377);
     rectangle(355,105,585,375);
     setcolor(14);
     if (ea mod 2)=0 then
     outtextxy(360,115,'       ОКНО ПОМОЩИ')
     else
     outtextxy(360,115,'       HELP WINDOW');
end;
procedure error1(ea:word);
begin
     setcolor(15);
     setfillstyle(1,12);
     bar(140,210,490,280);
     rectangle(145,215,485,275);
     rectangle(147,217,483,273);
     if (ea mod 2)=0 then
     begin
     outtextxy(150,227,'                 Ошибка!                 ');
     outtextxy(150,237,'       Вводимые параметр не число!!      ');
     outtextxy(150,250,' Проверьте значение и заново введите его.');
     end
     else
     begin
     outtextxy(150,227,'                Error!                   ');
     outtextxy(150,237,'  The value you entered isn`t a quantity!!');
     outtextxy(150,250,'       Check it and put it in again.     ');
     end;
     sound(600);
     delay(4000);
     nosound;
     readln;
     readln;
end;
procedure error(ea:word);
{Процедура ошибки}
begin
     setcolor(15);
     setfillstyle(1,12);
     bar(140,210,490,260);
     rectangle(145,215,485,255);
     rectangle(147,217,483,253);
     if (ea mod 2)=0 then
     begin
     outtextxy(150,227,'                 Ошибка!');
     outtextxy(150,237,'    Недостаток вводимых параметров!!');
     end
     else
     begin
     outtextxy(150,227,'                 Error!');
     outtextxy(150,237,'       Not all parameters are set!');
     end;
     sound(600);
     delay(4000);
     nosound;
     readln;
end;
procedure newsctext(ea:word);
{Текст для процедуры newsc}
begin
if ea mod 2 =0 then
      begin
      settextstyle(0,0,1);
     setcolor(15);
     outtextxy(400,440,'Язык - Русский. ');
     outtextxy(400,450,'Версия 1.0 Последнее издание');
     outtextxy(400,460,'й Все права защищены.');
     end
     else
     begin
         settextstyle(0,0,1);
         setcolor(15);
         outtextxy(400,440,'Language - English.');
         outtextxy(400,450,'Version 1.0 Final release.');
         outtextxy(400,460,'й All rights reserved.');
     end;
end;
procedure newsc(ea:word);
{Процедура обновления экрана}
begin
     cleardevice;
     setfillstyle(10,8);
     floodfill(1,1,15);
     setcolor(0);
     setfillstyle(1,7);
     bar(80,10,580,80);
     rectangle(82,12,578,78);
     rectangle(85,15,575,75);
     settextstyle(0,0,2);
     setcolor(10);
     if ea mod 2 =0 then
     begin
     settextstyle(0,0,2);
     outtextxy(90,20,'    Вычисление интеграл ');
     outtextxy(90,50,'   методом Ньютона-Котеса.');
     newsctext(ea);
     end
     else
     begin
     settextstyle(3,0,2);
     outtextxy(90,20,'            Calculeting of integral');
     outtextxy(90,47,'        using the Newton-Cotes method.');
     newsctext(ea);
     end;
     settextstyle(0,0,1);
end;
procedure winwin1;
{Окно процедуры win1}
begin
     setfillstyle(1,7);
     bar(160,110,460,380);
     setcolor(0);
     rectangle(162,113,457,377);
     rectangle(165,115,455,375);
end;
procedure win1(ea:word);
{Вводное окно}
begin
     settextstyle(0,0,1);
     setcolor(10);
     if (ea mod 2)=0 then
     begin
     outtextxy(168,135,'Министерство Высшего образования РФ);
     outtextxy(168,150,'Московский Государственный Институт');
     outtextxy(168,160,'         Электронной Техники       ');
     outtextxy(168,170,'      (Технический лниверситет)    ');
     outtextxy(168,180,'             Лицей №1557           ');
     outtextxy(168,210,'          КУРСОВАЯ РАБО'А          ');
     outtextxy(168,230,'      «Вычисление интеграла        ');
     outtextxy(168,245,'      метедом Ньютона-Котеса»      ');
     outtextxy(158,270,'       Написал: Коноплев А.А.      ');
     outtextxy(158,285,'  Руководитель: доцент Колдаев В.Д.');
     end
     else
     begin
     outtextxy(168,135,'    Department of High Education   ');
     outtextxy(168,150,'     Moscow State Institute of     ');
     outtextxy(168,160,'         Electronic Technics       ');
     outtextxy(168,170,'        (Technics University)      ');
     outtextxy(168,180,'            Lyceum №1557           ');
     outtextxy(168,210,'            COURSE WORK            ');
     outtextxy(168,230,'     «Calculation of integral      ');
     outtextxy(168,245,'      by Newton-Cotes method»      ');
     outtextxy(158,270,'       Author: Konoplev A.A.       ');
     outtextxy(158,285,'  Supervisor:senior lecturer       ');
     outtextxy(158,300,'                    Koldaev V.D.   ');
     end;
end;
procedure win2(ea:word;var k:word);
{Окно выбора  способа  подсчета }
var
c:char;
x:integer;
f:string;
begin
     setcolor(2);
     setfillstyle(1,5);
     bar(70,200,340,330);
     rectangle(75,205,335,325);
     rectangle(77,207,333,323);
     settextstyle(0,0,0);
     setfillstyle(1,15);
     bar(80,250,330,270);
     setfillstyle(1,5);
     bar(80,285,330,305);
     if ea mod 2 =0 then
     begin
     outtextxy(77,220,'Выбирете способ задания значений');
     outtextxy(75,230,'            функции.    ');
     outtextxy(70,255,'        По таблице(в ручную)');
     outtextxy(70,295,'        По расчетам(автом т.)');
     end
     else
     begin
     outtextxy(77,220,' Choose a method of putting in');
     outtextxy(75,230,'    the values of function.   ');
     outtextxy(70,255,'       By the table(by hand)');
     outtextxy(70,295,'      By calculations(automat.)');
     end;
     helpwin(ea);
     if ea mod 2 =0 then
     begin
     outtextxy(360,140,'В этом способе необходимо');
     outtextxy(360,155,'самостоятельно вводить');
     outtextxy(360,170,'значения функции.');
     end
     else
     begin
     outtextxy(360,140,'In this method you have');
     outtextxy(360,155,'to put in values of ');
     outtextxy(360,170,'function by yourself.');
     end;
     x:=0;
     repeat
     if keypressed then
     begin
          c:=readkey;
          if (c=#80) or (c=#72) then
             x:=x+1;
          setfillstyle(1,15);
          if (x mod 2)=0 then
             begin
                  bar(80,250,330,270);
                  setfillstyle(1,5);
                  bar(80,285,330,305);
                  helpwin(ea);
                  if ea mod 2 =0 then
                  begin
                       outtextxy(360,140,'В этом способе необходимо');
                       outtextxy(360,155,'самостоятельно вводить');
                       outtextxy(360,170,'значения функции.');
                  end
                  else
                  begin
                       outtextxy(360,140,'In this method you have');
                       outtextxy(360,155,'to put in values of ');
                       outtextxy(360,170,'function by yourself.');
                  end;
             end
          else
              begin
                   bar(80,285,330,305);
                   setfillstyle(1,5);
                   bar(80,250,330,270);
                   helpwin(ea);
                   if ea mod 2 =0 then
                   begin
                   outtextxy(360,140,'В этом способе компьютер');
                   outtextxy(360,155,'сам вычесляет значения');
                   outtextxy(360,170,'функции по вводимой функции.');
                   end
                   else
                   begin
                   outtextxy(360,140,'In this method PC will');
                   outtextxy(360,155,'automaticly count the value');
                   outtextxy(360,170,'of function by the function');
                   outtextxy(360,185,'you enter  ');
                   end;
              end;
     setcolor(2);
     if ea mod 2 =0 then
     begin
     outtextxy(70,255,'        По таблице(в ручную)');
     outtextxy(70,295,'        По расчетам(автом т.)');
     end
     else
     begin
     outtextxy(70,255,'       By the table(by hand)');
     outtextxy(70,295,'      By calculations(automat.)');
     end;
     end;
     until c=#13;
k:=x mod 2;
end;
procedure wwod1(ea:word;var y:array of double;var n:integer;var a,b:real);
{Окно ручного ввода функции}
var
i,p:integer;
s,f:string;
p1:real;
c:char;

begin
     wwodn(ea,n);
     if n=0 then
        wwodn(ea,n);
     newsc(ea);
     wwodab(ea,a,b);
     helpwin(ea);
     if ea mod 2 =0 then
     begin
     outtextxy(360,140,'В этом окне необходимо');
     outtextxy(360,155,'постепенно вводить');
     outtextxy(360,170,'значения функции.');
     outtextxy(360,185,'после каждого ввода');
     outtextxy(360,200,'определенного значения');
     outtextxy(360,215,'нажмите ENTER.');
     end
     else
     begin
     outtextxy(360,140,'In this window you have');
     outtextxy(360,155,'to gradually enter the');
     outtextxy(360,170,'values of functions.');
     outtextxy(360,185,'After each enter press');
     outtextxy(360,200,'ENTER key.');
     end;
     setfillstyle(1,9);
     bar(40,200,330,300);
     rectangle(45,205,325,295);
     rectangle(47,207,323,293);
     if ea mod 2 =0 then
     outtextxy(56,227,'Введите 0 -е значение финкции:')
     else
     outtextxy(56,227,' Enter  0 -th value of function:');
     for i:=0 to n do
     begin
         setfillstyle(1,0);
         bar(137,250,180,273);
         gotoxy(19,17);
         setfillstyle(1,9);
         read(p1);
         y[i]:=p1;
         bar(120,227,134,240);
         str(i+1,s);
         outtextxy(120,227,s);
         bar(310,220,320,250);
     end;

end;
procedure wwod2(ea:word;var ea1:word;var n:integer;var a,b:real;var
st:string);
{Окно 2 меню автомат. подсчета}
var
   i:integer;
   c,k:char;
   x:longint;
   f:string;
begin
     repeat
     x:=-600000;
     if keypressed then
        c:=readkey;
     c:='t';
     newsc(ea);
     setfillstyle(1,15);
     bar(70,120,342,330);
     setcolor(12);
     rectangle(75,125,337,325);
     rectangle(77,127,335,323);
     settextstyle(0,0,0);
     setfillstyle(1,11);
     bar(80,170,330,190);
     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;
     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');
          outtextxy(360,155,' in the number of units');
     end;
     repeat
     if keypressed then
     begin
          c:=readkey;
          case c of
               #80:
                   x:=x-1;
               #72:
                   x:=x+1;
          end;
          setfillstyle(1,11);
          case (abs(x) mod 5) of
               0:
                 begin
                      bar(80,170,330,190);
                      setfillstyle(1,15);
                      bar(80,200,330,220);
                      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 put');
                      outtextxy(360,155,'in the number of units.');
                      end;
                 end;
               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 (n<=0)or(a=b)or(st='') then
            error(ea);
          4:
            halt;
     end;
     until (n>0)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)<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)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-ков на коэффициенты}
     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]