заработок на кликах
Вы еще не зарегестрированы на Uchit.net? Зачем?
Login: Pass:

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

реферат: Кибернетика

Оцените работу
всего оценок0 общий балл0
Зарегистрируйтесь


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


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

(Технический Университет)


Лицей №1557










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



“Вычисление интеграла  методом

Ньютона-Котеса”




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

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





Москва, 2001г.
















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


























   

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

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






















Пусть некоторая функция f(x) задана в уздах интерполяции:

         (i=1,2,3…,n) на отрезке [а,b] таблицей значений:



X0=a
X1
X2

XN=b

Y0=f(x0)

Y1=f(x1)

Y2=f(x2)

YN=f(xN)


Требуется найти значение интеграла    .

Для начала составим интерполяционный многочлен Лагранджа:



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




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










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



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

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

Положим, что

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

Теперь рассмотрим несколько примеров.






Пример 1.

   Вычислить с помощью метода Ньютона-Котаса:                            , при 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.

Вычислить при помощи метода Ньютона-Котеса

     , взяв 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

  1. Находим коэффициенты Ньютона Котеса:

H0=H5=0.065972 ;H1=H4=0.260417 ;H2=H3=0.173611 ;

4)Подставим значения в формулу и получим:

  


    Рассмотрим частные случаи формулы Ньйтона-Котеса.

Пусть n=1 тогда

H0=H1=0.5 и конечная формула примет вид:

Тем самым в качестве частного случая нашей формулы мы получили формулу трапеций.

Взяв n=3, мы получим

. Частный случай формулы Ньютона Котеса формула Симпсона








     Теперь произведем анализ алгоритма и рассмотрим  основной принцип работы программы.

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

Основной проблемой вычисления коэффициентов является интеграл от произведения множителей. Для его расчета необходимо:


А) посчитать коэффициенты при раскрытии скобок при q

(процедура mnogoclen)

Б) домножить их на 1/n , где n степень при q (процедура koef)

В) подставить вместо q значение n (функция integral)


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


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



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





















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























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
















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














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













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










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





Программа состоит из 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)<b) and(k1<>0) then

          begin

          if (abs(240-k1)>2) then

          begin

          if k1<240 then

             k1:=k1+1

          else

              k1:=k1-1;

          if c>7 then

          setfillstyle(6,3)

          else

              setfillstyle(1,3);

          floodfill(320-round((240+x1)/10)+x,k1,15);

          end;

          end;

    end;

    end;

    str(x1,f2);

    outtextxy(1,450,f2);

    if (n mod 2)=0 then

    for i:=-9 to 10 do

    begin

         settextstyle(2,0,2);

         setcolor(14);

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

         begin

            str((i*24+round((240+x1)/10))/c:2:2,f);

            p:=247;

            outtextxy(310+i*24,p,f);

            str(-i*24/c:2:2,f);

            outtextxy(330,240+i*24,f);

         end;

    end;

    for i:=-9 to 10 do

     begin

          setcolor(15);

          if ((r mod 2)=1) and (i<>0) then

          begin

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

              line(320+i*24,20,320+i*24,460);

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

              line(80,240+i*24,560,240+i*24);

          end;

     end;

    setcolor(15);

    d:=readkey;

    case d of

          #75:

              begin

              x1:=x1-30;

              x2:=x2-30;

              end;

          #77:

          begin

               x1:=x1+30;

               x2:=x2+30;

          end;


          #80:

              if c>1 then

                   c:=c-1;

          #72:

              c:=c+1;

          #71:

              c:=24;

          #79:

              n:=n+1;

          #83:

              r:=r+1;

          #82:

              v:=v+1;

          #73:

              begin

                   c:=24;

                   n:=0;r:=0;v:=0;x1:=-240;x2:=240;

              end;

          #59:

              hwg(ea);

    end;


    until d=#13;

end;

end.


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

==========МОДУЛЬ UNIT==========

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

{$N+}

Unit k_unit;

{Модуль нахождения интеграл  от многочлена q(q-1)..(q-i+1)(q-i-1)..(q-n),}

{где n-точность интеграла ,i-номер коофициента.                          }

interface

procedure rasposn(f:string;x:real;var ec:word;var t:real);

procedure hkoef(n:integer;var h:array of double);

procedure funktia(n:integer;a,b:real;var y:array of double;c:real;f:string);

procedure koef(w:array of double;n:integer;var e:array of double);

procedure mnogochlen(n,i:integer;var c:array of double);

function facktorial(n:integer):double;

function integral(w:array of double;n:integer):double;

function mainint(n:integer;a,b:real;y:array of double):double;

implementation

procedure rasposn(f:string;x:real;var ec:word;var t:real);

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

var

k:word;

begin

     k:=pos('x',f);

     if k<>0 then

     begin   {Распознавание функции}

         ec:=1; {Код ошибки}

         t:=x;

         k:=pos('abs(x)',f);

         if k<>0 then t:=abs(x);

         k:=pos('sin(x)',f);

         if k<>0 then t:=sin(x);

         k:=pos('cos(x)',f);

         if k<>0 then t:=cos(x);

         k:=pos('arctg(x)',f);

         if k<>0 then t:=arctan(x);

         k:=pos('sqr(x)',f);

         if k<>0 then t:=x*x;

         k:=pos('exp(x)',f);

         if k<>0 then t:=exp(x);

         k:=pos('cos(x)*x',f);

         if k<>0 then t:=cos(x)*x;

         k:=pos('ln(x)',f);

         if k<>0 then

            begin

            if x>0 then t:=ln(x)

            else

                t:=0;

            end;

         k:=pos('sqrt(x)',f);

         if k<>0 then

            if x>=0 then t:=sqrt(x)

            else t:=0;

         k:=pos('arcctg(x)',f);

         if k<>0 then t:=pi/2-arctan(x);

         k:=pos('sin(x)/x',f);

         if k<>0 then if x<>0 then t:=sin(x)/x;

     end

     else

         ec:=0;

end;

procedure funktia(n:integer;a,b:real;var y:array of double;c:real;f:string);

{Процедур  подсчет  Y-ков и распознавания функции}

var 

t,h,x:real;

k,i:integer;

es:word;

begin

     h:=(b-a)/n;

     for i:=0 to n do

     begin

         x:=(a+h*i)/c;

         rasposn(f,x,es,t);

         y[i]:=t;

     end;

end;

procedure koef(w:array of double;n:integer;var e:array of double);

{Изменение коофициентов для интеграла}

var

t:integer;

begin

     for t:=1 to n do

         e[t]:=w[t]/(n-t+2);

end;

procedure mnogochlen(n,i:integer;var c:array of double);

{процедура нахождения коофициентов при Q^n(q в степени n )}

var

k,j:integer;

d:array[1..100] of double;

begin

     d[1]:=1;

     for j:=1 to n do    

     begin  {Вычисление коэффициентов при раскрытии q*(q-1)*(q-2)*..*(q-n)}

          d[j+1]:=d[j]*j*(-1);

          if j>1 then

             for k:=j downto 2 do

                 d[k]:=d[k]+d[k-1]*j*(-1);

     end;

     c[1]:=d[1]; {Деление многочлена на (q-i) по схеме Горнера}

     for j:=1 to n+1 do

         c[j]:=i*c[j-1]+d[j];

     koef(c,n,c); {Изменение коэффициентов при интегрировании}

end;

function facktorial(n:integer):double;

{функция нахождения факториала }

var

t:integer;

s:double;

begin

     s:=1;

     if n=0 then

        s:=1

     else

         for t:=1 to n do

             s:=s*t;

     facktorial:=s;

end;

function integral(w:array of double;n:integer):double;

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

var

t,p:integer;

s,c:double;

begin

     s:=0;p:=n;

     for t:=0 to p+1 do

         s:=s+w[t]*exp((p-t+2)*ln(p)); {Подсчет интеграла}

     integral:=s;

end;

procedure hkoef(n:integer;var h:array of double);

{Процедура подсчета коэф. Ньютона-Котеса}

var

p,j,d,c,i:integer;

kq:array[0..20] of double;

s:array[0..20] of double;

begin

     p:=n;

     if (p mod 2)=1 then {Вычисление половины от всех вычислений коэффициентов}

        d:=round((p-1)*0.5)

     else

         d:=round(0.5*p);

     for i:=0 to n do

     begin

         mnogochlen(p,i,kq);

         s[i]:=integral(kq,p); {Формирование массива из интегралов}

     end;

     for i:=0 to d do

         begin

              if ((p-i) mod 2) = 0 then

                 c:=1

              else

                  c:=(-1);

              h[i]:=(c*s[i])/(facktorial(i)*facktorial(p-i)*p);

              h[p-i]:=h[i];

         end;

end;

function mainint(n:integer;a,b:real;y:array of double):double;

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

var

sum:double;

p,i:integer;

kq,h:array[0..20] of double;

begin

     p:=n;

     hkoef(n,h);

     sum:=0;

     for i:=0 to p do

         sum:=sum+h[i]*y[i]; {Сумма произведений y-ков на коэффициенты}

     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.



















       Рассмотрим результаты тестовых испытаний для функций 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) желаемая формула для распознования.

  1.    Вся помощь по вводу и работе с пограммой выводится в окне помощи.






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


                        


  



В.И. Грызлов «Турбо Паскаль 7.0» Москва: ДМК 2000г.

Данилина «Численные методы» Москва: Высшая школа 1978г.

 

 
Дружить
Uchit.net в социальных сетях