Министерство Высшего Образования РФ.
Московский Институт Электронной Техники
(Технический Университет)
Лицей №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 ec0 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 ec0 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 ec0 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 (n0)and(ab)and(st"")and((abs(x) mod 5)=3);
end;
procedure win3(ea:word;n:integer;a,b:real;int:double;f:string;h:array of
double;var k:word);
{Последнее окно просмотра результатов}
var i:integer; c:char; x:longint; p1,p:string; y:array[0..16] of double;
begin funktia(n,a,b,y,1,f); f:="("+f+")"+"dx ="; repeat x:=-600000; newsc(ea); setfillstyle(1,2); bar(170,120,490,360); setcolor(14); rectangle(175,125,485,355); rectangle(177,127,483,353); settextstyle(0,0,0); setfillstyle(1,1); bar(180,170,480,190); if ea mod 2 =0 then begin outtextxy(180,135,Функция распознана.Интеграл подсчитан."); outtextxy(180,180," Посмотреть значение интеграла"); outtextxy(180,210,"Посмотреть коэффициенты Ньютона-Котеса"); outtextxy(180,240," Посмотреть значения функции"); outtextxy(180,270," Посмотреть график" ); outtextxy(180,300," Считать снова"); outtextxy(180,330," Выход "); end else begin outtextxy(180,135,"Function Indentified.Integral counted."); outtextxy(180,180," View value of integral"); outtextxy(180,210," View Newton-Cotes coefficients"); outtextxy(180,240," Veiw values of function"); outtextxy(180,270," View graphik " ); outtextxy(180,300," Count again"); outtextxy(180,330," Exit "); end; repeat if keypressed then begin c:=readkey; case c of
#80: x:=x-1;
#72: x:=x+1; end; setfillstyle(1,1); case (abs(x) mod 6) of
0: begin bar(180,170,480,190); setfillstyle(1,2); bar(180,200,480,220); bar(180,320,480,340); end;
1: begin bar(180,200,480,220); setfillstyle(1,2); bar(180,170,480,190); bar(180,230,480,250); end;
2: begin bar(180,230,480,250); setfillstyle(1,2); bar(180,200,480,220); bar(180,260,480,280); end;
3: begin bar(180,260,480,280); setfillstyle(1,2); bar(180,230,480,250); bar(180,290,480,310); end;
4: begin bar(180,290,480,310); setfillstyle(1,2); bar(180,260,480,280); bar(180,320,480,340); end;
5: begin bar(180,320,480,340); setfillstyle(1,2); bar(180,290,480,310); bar(180,170,480,190); end; end; if ea mod 2 =0 then begin outtextxy(180,135,"Функция распознана.Интеграл подсчитан."); outtextxy(180,180," Посмотреть значение интеграла"); outtextxy(180,210,"Посмотреть коэффициенты Ньютона-Котеса"); outtextxy(180,240," Посмотреть значения функции"); outtextxy(180,270," Посмотреть график " ); outtextxy(180,300," Считать снова"); outtextxy(180,330," Выход "); end else begin outtextxy(180,135,"Function Indentified.Integral counted."); outtextxy(180,180," View value of integral"); outtextxy(180,210," View Newton-Cotes coefficients"); outtextxy(180,240," Veiw values of function"); outtextxy(180,270," View graphik " ); outtextxy(180,300," Count again"); outtextxy(180,330," Exit "); end;
end; until c=#13; c:="t"; case (abs(x) mod 6) of
0:begin setcolor(15); setfillstyle(1,12); bar(140,200,490,280); rectangle(145,205,485,275); rectangle(147,207,483,273); settextstyle(2,0,1); setusercharsize(1,1,5,1); outtextxy(170,210,"S"); settextstyle(2,0,4); str(a:3:3,p); outtextxy(160,257,p); str(b:3:3,p); outtextxy(160,212,p); settextstyle(3,0,2); outtextxy(180,224,f); p:=""; str(abs(int):7:3,p); outtextxy(190+length(f)*12,224,p); readln; end;
1: begin newsc(ea); setfillstyle(1,2); bar(170,120,490,180+n*15); setcolor(14); rectangle(175,125,485,175+n*15); rectangle(177,127,483,173+n*15); if ea mod 2 =0 then begin outtextxy(180,130,"Коэффициенты Ньютона-Котеса:"); outtextxy(180,140+(n+1)*15,"Нажмите ENTER для продолжения"); end else begin outtextxy(180,130,"Newton-Cotes coefficients:"); outtextxy(180,140+(n+1)*15,"Press ENTER to continue"); end; hkoef(n,h); for i:=0 to n do begin str(i,p);str(h[i]:2:4,p1); p:="H"+p+" = "+p1; outtextxy(180,140+i*15,p); end; readln; end;
2:begin newsc(ea); setfillstyle(1,2); bar(170,120,490,180+n*15); setcolor(14); rectangle(175,125,485,175+n*15); rectangle(177,127,483,173+n*15); if ea mod 2 =0 then begin outtextxy(180,130,"Значения функции:"); outtextxy(180,140+(n+1)*15,"Нажмите ENTER для продолжения"); end else begin outtextxy(180,130,"Values of function:"); outtextxy(180,140+(n+1)*15,"Press ENTER to continue"); end; for i:=0 to n do begin str(i,p);str(y[i]:2:4,p1); p:="Y"+p+" = "+p1; p1:=""; outtextxy(180,140+i*15,p); str((a+i*(b-a)/n):2:4,p1); str(i,p); if ea mod 2 = 0 then p:=",При "+"X"+p+" = "+p1 else p:=",When "+"X"+p+" = "+p1; outtextxy(285,140+i*15,p); end;
readln; end;
3: graphik(ea,a,b,f);
5: begin closegraph; halt; end; end; until (abs(x) mod 6)=4; k:=abs(x) mod 6; end; end.
================================================
========МОДУЛЬ GRAPHIC========
================================================
unit graphic;
interface
uses
k_unit,crt,graph;
procedure hwg(ea:word);
procedure graphik(ea:word;a,b:real;f1:string);
implementation
procedure hwg(ea:word);
{Процедура окна помощи при графике}
var
f:string;
begin settextstyle(0,0,0); setfillstyle(1,3); bar(150,100,390,380); setcolor(0); rectangle(153,103,387,377); rectangle(155,105,385,375); setcolor(14); if ea mod 2 =0 then begin outtextxy(160,115," ОКНО ПОМОЩИ"); outtextxy(160,140," Для работы с графиком"); outtextxy(160,155," используйте клавиши:"); outtextxy(160,180," PAGE UP-первоначальный"); outtextxy(160,195," вид графика;"); outtextxy(160,210," HOME-начальный масштаб;"); outtextxy(160,225," INSERT-включить/выключеть"); outtextxy(160,240," заливку области;"); outtextxy(160,255," DELETE-включить/выключеть"); outtextxy(160,270," сетку;"); outtextxy(160,285," END-показать/убрать цифры"); outtextxy(160,300," F1- Помощь;"); outtextxy(160,315," Стрелки ВВЕРХ/ВНИЗ- "); outtextxy(160,330," увеличение/уменьшение"); outtextxy(160,345," масштаб ."); outtextxy(160,360,"Для возрата нажмите ENTER."); end else begin outtextxy(160,115," HELP WINDOW"); outtextxy(160,140," For the work with graphic"); outtextxy(160,155," use this keys:"); outtextxy(160,180," PAGE UP-Primery form of"); outtextxy(160,195," graphik;"); outtextxy(160,210," HOME-Primery scale;"); outtextxy(160,225," INSERT-Turn on/off inking"); outtextxy(160,240," the field;"); outtextxy(160,255," DELETE-Turn on/off the"); outtextxy(160,270," net;"); outtextxy(160,285," END-View/delete the figures"); outtextxy(160,300," F1- Help;"); outtextxy(160,315," Arrows UP/DOWN-Increase/ "); outtextxy(160,330," lower the scale;"); outtextxy(160,360,"Press ENTER to continue."); end; readln; setcolor(15);
end;
procedure graphik(ea:word;a,b:real;f1:string);
{процедура построения графиков}
var
f,f2:string;
d:char;
i,v,r:integer;
x1,x2,n,p,x:integer;
c,k,k1:longint;
y:array[0..1] of double;
begin x1:=-240; x2:=240; c:=24; setcolor(15); n:=0;v:=0;r:=0; repeat cleardevice; settextstyle(0,0,0); if ea mod 2 =0 then begin outtextxy(10,1,"Нажмите F1 для помощи"); str(c/24:2:2,f); f:="Масштаб "+f+":1"; end else begin outtextxy(10,1,"Press F1 for help"); str(c/24:2:2,f); f:="Scale "+f+":1"; end; outtextxy(200,1,f); settextstyle(3,0,1); outtextxy(307,10,"y"); outtextxy(574,235,"x"); outtextxy(310,240,"0"); setlinestyle(1,7,100); line(70,240,580,240); line(320,20,320,460); line(320,20,315,25); line(321,20,326,25); line(580,239,575,244); line(580,240,575,235); line(70,239,580,239); line(321,20,321,460); for i:=-9 to 10 do begin if ((320+i*24)71) then line(320+i*24,240,320+i*24,242); if ((240+i*24)19) then line(320,240+i*24,322,240+i*24); end; setcolor(15); for x:= -240+round((240+x1)/10) to 240+round((240+x1)/10) do begin funktia(1,x-1,x,y,c,f1); k:=round(240-(y[0])*c); k1:=round(240-(y[1])*c); if ((k0)or(k10)) then line(319-round((240+x1)/10)+x,k,320-round((240+x1)/10)+x,k1); end; if (v mod 2)=0 then begin funktia(1,a,b,y,1,f1); k:=round(240-(y[0])*c); k1:=round(240-(y[1])*c); line(320-round((240+x1)/10)+round(a*c),k,320-
round((240+x1)/10)+round(a*c),240); line(320-round((240+x1)/10)+round(b*c),k1,320-
round((240+x1)/10)+round(b*c),240); if 320-round((240+x1)/10)+a*c560 then begin funktia(1,(-240-round((240+x1)/10))/c,(240-
round((240+x1)/10))/c,y,1,f1); k1:=round(240-(y[1])*c); line(560,k1,560,240); end; for x:= -240 to 240 do begin funktia(1,x-1,x,y,c,f1); k1:=round(240-(y[1])*c); if ((x/c)>a) and ((x/c)2) then begin if k17 then setfillstyle(6,3) else setfillstyle(1,3); floodfill(320-round((240+x1)/10)+x,k1,15); end; end; end; end; str(x1,f2); outtextxy(1,450,f2); if (n mod 2)=0 then for i:=-9 to 10 do begin settextstyle(2,0,2); setcolor(14); if ((320+i*24)71)and(i0) then begin str((i*24+round((240+x1)/10))/c:2:2,f); p:=247; outtextxy(310+i*24,p,f); str(-i*24/c:2:2,f); outtextxy(330,240+i*24,f); end; end; for i:=-9 to 10 do begin setcolor(15); if ((r mod 2)=1) and (i0) then begin if ((320+i*24)71) then line(320+i*24,20,320+i*24,460); if ((240+i*24)19) then line(80,240+i*24,560,240+i*24); end; end; setcolor(15); d:=readkey; case d of
#75: begin x1:=x1-30; x2:=x2-30; end;
#77: begin x1:=x1+30; x2:=x2+30; end;
#80: if c>1 then c:=c-1;
#72: c:=c+1;
#71: c:=24;
#79: n:=n+1;
#83: r:=r+1;
#82: v:=v+1;
#73: begin c:=24; n:=0;r:=0;v:=0;x1:=-240;x2:=240; end;
#59: hwg(ea); end;
until d=#13; end; end.
================================================
==========МОДУЛЬ UNIT==========
================================================
{$N+}
Unit k_unit;
{Модуль нахождения интеграл от многочлена q(q-1)..(q-i+1)(q-i-1)..(q-n),}
{где n-точность интеграла ,i-номер коофициента. }
interface
procedure rasposn(f:string;x:real;var ec:word;var t:real);
procedure hkoef(n:integer;var h:array of double);
procedure funktia(n:integer;a,b:real;var y:array of
double;c:real;f:string);
procedure koef(w:array of double;n:integer;var e:array of double);
procedure mnogochlen(n,i:integer;var c:array of double);
function facktorial(n:integer):double;
function integral(w:array of double;n:integer):double;
function mainint(n:integer;a,b:real;y:array of double):double;
implementation
procedure rasposn(f:string;x:real;var ec:word;var t:real);
{Процедура распознования функции}
var
k:word;
begin k:=pos("x",f); if k0 then begin {Распознавание функции} ec:=1; {Код ошибки} t:=x; k:=pos("abs(x)",f); if k0 then t:=abs(x); k:=pos("sin(x)",f); if k0 then t:=sin(x); k:=pos("cos(x)",f); if k0 then t:=cos(x); k:=pos("arctg(x)",f); if k0 then t:=arctan(x); k:=pos("sqr(x)",f); if k0 then t:=x*x; k:=pos("exp(x)",f); if k0 then t:=exp(x); k:=pos("cos(x)*x",f); if k0 then t:=cos(x)*x; k:=pos("ln(x)",f); if k0 then begin if x>0 then t:=ln(x) else t:=0; end; k:=pos("sqrt(x)",f); if k0 then if x>=0 then t:=sqrt(x) else t:=0; k:=pos("arcctg(x)",f); if k0 then t:=pi/2-arctan(x); k:=pos("sin(x)/x",f); if k0 then if x0 then t:=sin(x)/x; end else ec:=0;
end;
procedure funktia(n:integer;a,b:real;var y:array of
double;c:real;f:string);
{Процедур подсчет Y-ков и распознавания функции}
var
t,h,x:real;
k,i:integer;
es:word;
begin h:=(b-a)/n; for i:=0 to n do begin x:=(a+h*i)/c; rasposn(f,x,es,t); y[i]:=t; end;
end;
procedure koef(w:array of double;n:integer;var e:array of double);
{Изменение коофициентов для интеграла}
var
t:integer;
begin for t:=1 to n do e[t]:=w[t]/(n-t+2);
end;
procedure mnogochlen(n,i:integer;var c:array of double);
{процедура нахождения коофициентов при Q^n(q в степени n )}
var
k,j:integer;
d:array[1..100] of double;
begin d[1]:=1; for j:=1 to n do begin {Вычисление коэффициентов при раскрытии q*(q-1)*(q-2)*..*(q-n)} d[j+1]:=d[j]*j*(-1); if j>1 then for k:=j downto 2 do d[k]:=d[k]+d[k-1]*j*(-1); end; c[1]:=d[1]; {Деление многочлена на (q-i) по схеме Горнера} for j:=1 to n+1 do c[j]:=i*c[j-1]+d[j]; koef(c,n,c); {Изменение коэффициентов при интегрировании}
end;
function facktorial(n:integer):double;
{функция нахождения факториала }
var
t:integer;
s:double;
begin s:=1; if n=0 then s:=1 else for t:=1 to n do s:=s*t; facktorial:=s;
end;
function integral(w:array of double;n:integer):double;
{функция подсчета самого интеграла}
var
t,p:integer;
s,c:double;
begin s:=0;p:=n; for t:=0 to p+1 do s:=s+w[t]*exp((p-t+2)*ln(p)); {Подсчет интеграла} integral:=s;
end;
procedure hkoef(n:integer;var h:array of double);
{Процедура подсчета коэф. Ньютона-Котеса}
var
p,j,d,c,i:integer;
kq:array[0..20] of double;
s:array[0..20] of double;
begin p:=n; if (p mod 2)=1 then {Вычисление половины от всех вычислений
коэффициентов} d:=round((p-1)*0.5) else d:=round(0.5*p); for i:=0 to n do begin mnogochlen(p,i,kq); s[i]:=integral(kq,p); {Формирование массива из интегралов} end; for i:=0 to d do begin if ((p-i) mod 2) = 0 then c:=1 else c:=(-1); h[i]:=(c*s[i])/(facktorial(i)*facktorial(p-i)*p); h[p-i]:=h[i]; end;
end;
function mainint(n:integer;a,b:real;y:array of double):double;
{функция подсчета основного интеграла}
var
sum:double;
p,i:integer;
kq,h:array[0..20] of double;
begin p:=n; hkoef(n,h); sum:=0; for i:=0 to p do sum:=sum+h[i]*y[i]; {Сумма произведений y-ков на коэффициенты} mainint:=sum*(b-a);
end;
end.
================================================
=======ОСНОВНАЯ ПРОГРАММА=======
================================================
{$N+} program Newton_Cotes_metod;{Программа нахождения определенного интеграла} uses {методом Ньютона-Котеса } k_unit,k_graph,graph,crt; const t=15; var c:char; a1,b1,a,b:real; n1,v,r,n:integer; h,y:array[0..t] of double; ea,k:word; int:double; f:string; begin ea:=10; v:=detect; initgraph(v,r,""); cleardevice; newsc(ea); winwin1; setcolor(15); outtextxy(380,430,"Нажмите F2 для смены языка."); repeat win1(ea); settextstyle(3,0,1); outtextxy(178,340,"Press Enter..."); delay(13000); bar(178,340,350,365); delay(13000); if keypressed then {Смена языка} begin c:=readkey; if c=#60 then begin ea:=ea+1; newsc(ea); winwin1; setcolor(15); if ea mod 2 =0 then outtextxy(380,430,"Нажмите F2 для смены языка.") else outtextxy(380,430,"Press F2 key to change language."); end; end; until c=#13; repeat newsc(ea); win2(ea,k); {Ввод способа задания функции} case k of
0: wwod1(ea,y,n,a,b);
1: begin wwod2(ea,ea,n1,a1,b1,f); n:=n1;a:=a1;b:=b1; k:=4; end; end; if k=4 then funktia(n,a,b,y,1,f); int:=mainint(n,a,b,y); {Вычисление интеграла} hkoef(n,h); proline(ea); win3(ea,n,a,b,int,f,h,k); {Последнее меню вывода результатов} until k4; closegraph; end.
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
Рассмотрим результаты тестовых испытаний для функций sin(x) на интервале [-5;3] и exp(x) на интервале [2;8]
| |n=1 |n=2 |n=3 |n=4 |n=5 |n=7 |
|Sin(x) |4,040017 |3,02112 |0,087629 |1,779012 |1,537481 |1,246 |
|Exp(x) |8965,041 |3581,999 |3271,82 |3002,908 |2990,644 |2974,322 |
|N=9 |n=12 |
|1,273561 |1,27366 |
|2973,593 |2973,569 |
Видно, что при увеличении числа узлов интерполяции точность растет,
однако при больших n (n>15) наблюдался обратный эффект.
Рекомендуемый диапозон n: от 7 до 13.
1) Интерфейс программы составлен на 2 языках: русском и английском. Переход с одного языка на другой осуществляется в вводном окне путем нажатия клавишы F2. Сменить язык можно только в этой части программы.
2) При вводе значений функции вручную необходимо вводить только цифры и после каждого ввода нажимать клавишу ENTER.
3) При испытании программы под разные операционные системы(Dos, Windows 98-
2k,NT, из под паскаля) происходил непонятный баг с неверным выводом на экран значений коэффициентов Ньютона-Котеса, хотя интеграл считался верно. Для нормального нахождения их желательно запускать программу через
Dos.
4) При вводе параметров в “Меню задания параметров нахождения интеграла” желательно их вводить постепенно сверху вниз, т.е. сначала ввести количество узлов интерполяции, затем пределы интегрирования, а уж потом вводить саму функцию.
5) Данная версия программы не способна распознавать все функции. Она может распознать только стандартные функции Турбо Паскаля и еще несколько дотполнительных: sin(x)/x, cos(x)*x ,arcctg(x). Для работы со специфическими функциями необходимо в модуле K-unit в процедуре RASPOSN в конце, перед end else, добавить : k:=pos(‘Формула f(x)^,f); if k0 then t:= ‘Формула f(x)^;
где ‘Формула f(x)^ – желаемая формула для распознования.
6) Вся помощь по вводу и работе с пограммой выводится в окне помощи.
Для нахождения интеграла существует много методов, однако, метод
Ньютона-Котеса один из самых быстрых: достаточно знать значения
коэффициентов для n=4, чтобы с точностью до сотых мгновенно посчитать
интеграл. Быстрота и простота –главные части этого метода.
В.И. Грызлов «Турбо Паскаль 7.0» Москва: ДМК 2000г.
Данилина «Численные методы» Москва: Высшая школа 1978г.
-----------------------
[pic]
[pic]
(zip - application/zip)
Статистика файлового архива
