Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)
Министерство Высшего Образования РФ.
Московский Институт Электронной Техники
(Технический Университет)
Лицей №1557
КУРСОВАЯ РАБОТА
“Вычисление интеграла методом
Ньютона-Котеса”
Написал: Коноплев А.А.
Проверил: доцент Колдаев В.Д.
Москва, 2001г.
Введение..................................................................................... 3
Теоретическая часть...................................................................4
Алгоритм работы........................................................................8
Код программы.........................................................................17
Модуль K_graph............................................................17
Модуль Graphic.............................................................34
Модуль K_unit...............................................................38
Основная программа....................................................40
Тестовые испытания.................................................................42
Полезные советы по работе с программой.............................42
Окна ввода и вывода программы.............................................
Вывод..........................................................................................43
Список литературы...................................................................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;
Вычисление:
Определим шаг h=(8-4)/5=0.8
Найдем значения 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 |
Находим коэффициенты Ньютона –Котеса:
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 (не обязательно) – файл установки русского языка.
Для работы программы с русским интерфайсом желательно запускать ее в режиме 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;