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

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

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

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

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

Лицей №1557

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

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

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

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

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

Москва, 2001г.

1.

Введение..................................................................

................... 3

2. Теоретическая

часть...................................................................4

3. Алгоритм

работы....................................................................

....8

4. Код

программы.................................................................

........17

. Модуль

K_graph............................................................17

. Модуль

Graphic.............................................................34

. Модуль

K_unit...............................................................38

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

5. Тестовые

испытания.................................................................

42

6. Полезные советы по работе с программой.............................42

7. Окна ввода и вывода

программы.............................................

8.

Вывод.....................................................................

.....................43

9. Список

литературы................................................................

...44

Математика - одна из самых древних наук. Труды многих ученых вошли в

мировой фонд и стали основой современных алгебры и геометрии. В конце XVII

в., когда развитие науки шло быстрыми темпами, появились понятия

дифференцирование, а вслед за ним и интегрирование. Многие правила

нахождения неопределенного интеграла в то время не были известны, поэтому

ученые пытались найти другие, обходные пути поиска значений. Первым

методом явился метод Ньютона – поиск интеграла через график функции, т.е.

нахождение площади под графиком, методом прямоугольников, в последствии

усовершенствованный в метод трапеций. Позже был придуман параболический

метод или метод Симпсона. Однако часть ученых терзал вопрос: А можно ли

объединить все эти методы в один??

Ответ на него был дан одновременно двумя математиками Ньютоном и Котесом.

Они вывели общую формулу, названную в их честь. Однако их метод был

частично забыт. В этой работе будут изложены основные положения теории,

рассмотрены различные примеры, приведены таблицы, полученные при различных

погрешностях, и конечно описана работа и код программы, рассчитывающей

интеграл методом Ньютона-Котеса.

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

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

|X0=a |X1 |X2 |… |XN=b |

|Y0=f(x0) |Y1=f(x1) |Y2=f(x2) |… |YN=f(xN) |

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

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

[pic]

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

[pic]

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

интерполяционным многочленом Лагранжа:

[pic]

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

постоянные элементы:

[pic]

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

[pic]

Для равноотстоящих узлов интерполяции на отрезке [a,b] величина шаг

определяется как h=(a-b)/n. Представив это выражение для h в формулу (4) и

вынося (b-a) за знак суммы, получим:

[pic]

Положим, что

[pic]

где i=0,1,2…,n; Числа Hi называют коэффициентами Ньютона-Котеса. Эти

коэффиценты не зависят от вида f(x), а являются функцией только по n.

Поэтому их можно вычислить заранее. Окончательная формула выглядит так:

[pic]

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

Пример 1.

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

, при n=7.

Вычисление.

1) Определим шаг: h=(7-0)/7=1.

2)Найдем значения y:

|x0=0 |y0=1 |

|x1=1 |y1=0.5 |

|x2=2 |y2=0.2 |

|x3=3 |y3=0.1 |

|x4=4 |y4=0.0588 |

|x5=5 |y5=0.0384 |

|x6=6 |y6=0.0270 |

|x7=7 |y7=0.02 |

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

H1=H7=0.0435, H1=H6=0.2040, H2=H5=0.0760 ,H3=H4=0.1730

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

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

Пример 2.

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

[pic] , взяв n=5;

Вычисление:

1) Определим шаг h=(8-4)/5=0.8

2) Найдем значения y:

|x0=0 |y0=-2.61 |

|x1=4.8 |y1=0.42 |

|x2=5.6 |y2=4.34 |

|x3=6.4 |y3=6.35 |

|x4=7.2 |y4=4.38 |

|x5=8 |y5=-0.16 |

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

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

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

[pic]

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

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

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

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

формулу трапеций.

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

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

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

работы программы.

Для вычисления интеграла сначала находятся коэффициенты Ньютона-

Котеса. Их нахождение осуществляется в процедуре hkoef.

Основной проблемой вычисления коэффициентов является интеграл от

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

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

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

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

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

Далее вычисляем факториалы (функция faktorial) и перемножаем

полученные выражения (функция mainint). Для увеличения быстроты работы

вводится вычисление половины от количества узлов интерполяции и последующей

подстановкой их вместо неподсчитанных.

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

[pic]

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

[pic]

[pic]

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

[pic]

[pic]

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

[pic]

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

[pic]

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

[pic]

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

[pic]

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

[pic]

Программа состоит из 8 файлов:

. K_main.exe – файл загрузки основной программы

. K_unit.tpu – модуль вычислительных процедур и функций

. K_graph.tpu – модуль графических процедур

. Graphic.tpu – модуль процедур для построения графика

. Egavga.bgi – файл графической инициализации

. Sans.chr, litt.chr – файлы шрифтов

. Keyrus.com (не обязательно) – файл установки русского языка.

Для работы программы с русским интерфайсом желательно запускать ее в

режиме DOS.

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

==========МОДУЛЬ GRAPH==========

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

{$N+}

unit k_graph;

interface

uses

crt,graph,k_unit,graphic;

procedure winwin1;

procedure proline(ea:word);

procedure winwwodab(ea:word);

procedure error1(ea:word);

procedure helpwin(ea:word);

procedure error(ea:word);

procedure newsctext(ea:word);

procedure newsc(ea:word);

procedure win1(ea:word);

procedure win2(ea:word;var k:word);

procedure wwodn(ea:word;var n:integer);

procedure wwodab(ea:word;var a,b:real);

procedure wwod1(ea:word;var y:array of double;var n:integer;var a,b:real);

procedure wwod2(ea:word;var ea1:word;var n:integer;var a,b:real;var

st:string);

procedure win3(ea:word;n:integer;a,b:real;int:double;f:string;h:array of

double;var k:word);

implementation

procedure proline(ea:word);

{Проседура полосы процесса}

var

i:integer;

f:string;

c:char;

begin

newsc(ea);

setcolor(15);

setfillstyle(1,7);

bar(160,150,460,260);

rectangle(165,155,455,255);

rectangle(167,157,453,253);

case (ea mod 2) of

0: outtextxy(180,170,' Идет работа .Ждите..');

1: outtextxy(180,170,' Working.Please wait..');

end;

setfillstyle(1,12);

setcolor(0);

rectangle(200,199,401,221);

for i:=1 to 9 do

line(200+i*20,200,200+i*20,220);

delay(20000);

for i:=1 to 100 do

begin

if ((i-1) mod 10)=0 then

line(200+((i-1) div 10)*20,200,200+((i-1) div 10)*20,220);

bar(round(200+2*(i-0.5)),200,200+2*i,220);

delay(1100);

setcolor(15);

setfillstyle(1,7);

bar(280,230,323,250);

str(i,f);

f:=f+'%';

outtextxy(290,235,f);

if (i mod 25) =0 then

bar(170,180,452,198);

if (ea mod 2)=0 then

case (i div 25) of

0:

outtextxy(170,190,'Подготовка ');

1:

outtextxy(170,190,'Расчет коеффициентов в многочлене');

2:

outtextxy(170,190,'Расчет коеффициентов Ньютона-Котеса');

3:

outtextxy(170,190,'Расчет интеграла');

end

else

case (i div 25) of

0:

outtextxy(170,190,'Prepearing');

1:

outtextxy(170,190,'Calculation of mnogochlen coeff.');

2:

outtextxy(170,190,'Calculation of Newton-Cotes coeff. ');

3:

outtextxy(170,190,'Calculation of integral');

end;

setfillstyle(1,12);

setcolor(0);

end;

end;

procedure winwwodn(ea:word);

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

var

c:char;

f:string;

begin

helpwin(ea);

if (ea mod 2) =0 then

begin

outtextxy(360,140,' В этом окне необходимо ');

outtextxy(360,155,' ввести количество узлов ');

outtextxy(360,170,' интерполяции, от которого ');

outtextxy(360,185,' будет зависить точность ');

outtextxy(360,200,' вычисления интеграл и ');

outtextxy(360,215,' количество зн чений функции.');

outtextxy(360,240,' ВНИМАНИЕ : НАСТОЯТЕЛЬНО ');

outtextxy(360,250,' РЕКОМЕНДУЕТСЯ НЕ ВВОДИТЬ ');

outtextxy(360,260,' ЗНАЧЕНИЕ N БОЛЬШЕ 12 !! ');

end

else

begin

outtextxy(360,140,' In this window you have to ');

outtextxy(360,155,' put into the number. ');

outtextxy(360,170,' The accuracy of calculation ');

outtextxy(360,185,' and the number of function ');

outtextxy(360,200,' parameters will depend on ');

outtextxy(360,215,' this number. ');

outtextxy(360,240,' WARNING: IT IS HARDLY ');

outtextxy(360,250,' RECOMENDED NOT TO PUT IN ');

outtextxy(360,260,' NUMBER MORE THEN 12 !! ');

end;

setcolor(2);

setfillstyle(1,14);

bar(70,200,340,300);

rectangle(75,205,335,295);

rectangle(77,207,333,293);

if (ea mod 2) =0 then

begin

outtextxy(90,227,'Введите количество узлов(n):');

outtextxy(80,270,'ВНИМАНИЕ: При больших n возможна');

outtextxy(80,280,'некорректная работа компьютера!!');

end

else

begin

outtextxy(80,217,'Put in number of');

outtextxy(80,227,' interpolation units:');

outtextxy(80,270,'WARNING:if you use big number ');

outtextxy(80,280,'of units,PC wont work properly!');

end;

setfillstyle(1,0);

bar(190,240,230,255);

end;

procedure wwodn(ea:word;var n:integer);

{Процедура ввода узлов n}

var

ec,p:integer;

k,f:string;

x:integer;

c:char;

begin

newsc(ea);

winwwodn(ea);

repeat

repeat

winwwodn(ea);

gotoxy(25,16);

read(k);

val(k,p,ec);

if ec<>0 then

begin

error1(ea);

readln;

end;

until ec=0;

n:=p;

if n>12 then

begin

if keypressed then

c:=readkey;

c:='r';

setcolor(15);

setfillstyle(1,12);

bar(140,210,490,300);

rectangle(145,215,485,295);

rectangle(147,217,483,293);

if (ea mod 2) =0 then

begin

outtextxy(150,227,' Предупреждение!');

outtextxy(150,237,' Вы дейcтвительно хотите использовать');

outtextxy(150,250,' большое значение N ???');

end

else

begin

outtextxy(150,227,' Warning!!

');

outtextxy(150,237,' Do you realy want to use a big

');

outtextxy(150,250,' number interpolation units(N)???

');

end;

sound(600);

delay(4000);

nosound;

setfillstyle(1,2);

bar(320,260,350,280);

setfillstyle(1,12);

bar(250,260,280,280);

repeat

if keypressed then

begin

c:=readkey;

if (c=#80) or (c=#72) or (c=#77) or (c=#75) then

x:=x+1;

setfillstyle(1,2);

if (x mod 2)=0 then

begin

bar(250,260,280,280);

setfillstyle(1,12);

bar(320,260,350,280);

end

else

begin

bar(320,260,350,280);

setfillstyle(1,12);

bar(250,260,280,280);

END;

end;

if (ea mod 2) =0 then

begin

outtextxy(255,267,'ДА');

outtextxy(325,267,'НЕТ');

end

else

begin

outtextxy(255,267,'YES');

outtextxy(325,267,'NO');

end;

until c=#13;

if abs(x mod 2)=1 then

begin

n:=0;

setcolor(15);

setfillstyle(1,2);

bar(160,200,460,280);

rectangle(165,205,455,275);

rectangle(167,207,453,273);

if (ea mod 2)=0 then

begin

outtextxy(180,227,'Для работы программы необходимо');

outtextxy(180,237,' заново ввести N.');

outtextxy(180,247,' Нажмите ENTER для продолжения.');

end

else

begin

outtextxy(180,227,' To continue you have to ');

outtextxy(180,237,' again put in N. ');

outtextxy(180,247,' Press ENTER to continue.');

end;

readln;

readln;

end;

end;

until n>0;

end;

procedure winwwodab(ea:word);

{Окно ввода приделов интегрирования}

var

f:string;

begin

helpwin(ea);

if (ea mod 2)=0 then

begin

outtextxy(360,140,' В этом окне необходимо');

outtextxy(360,155,' ввести сначала нижнее');

outtextxy(360,170,' значение интеграл и нажать');

outtextxy(360,185,' ENTER, а затем ввести');

outtextxy(360,200,' верхнее значение интеграла');

outtextxy(360,215,' и снова нажать ENTER.');

end

else

begin

outtextxy(360,140,' In this window you have to:');

outtextxy(360,155,'firstly, put in lower value ');

outtextxy(360,170,'of integral and press ENTER,');

outtextxy(360,185,'then put in higher value');

outtextxy(360,200,'of integral and press ENTER');

end;

setcolor(2);

setfillstyle(1,5);

bar(10,210,335,320);

rectangle(15,215,330,315);

rectangle(17,217,328,313);

settextstyle(0,0,0);

if (ea mod 2)=0 then

begin

outtextxy(20,230,' Введите нижнее значение');

outtextxy(20,244,' интеграл :');

outtextxy(20,262,' Введите верхнее значение');

outtextxy(20,272,'интеграл :');

end

else

begin

outtextxy(20,230,' Put in lower value of');

outtextxy(20,244,' integral:');

outtextxy(20,262,' Put in higher value of');

outtextxy(20,272,'integral:');

end;

end;

procedure wwodab(ea:word;var a,b:real);

{Процедура ввода приделов интегрирования}

var

f:string;

k:string;

ec:integer;

begin

newsc(ea);

winwwodab(ea);

readln;

repeat

winwwodab(ea);

gotoxy(16,16);

read(k);

val(k,a,ec);

if ec<>0 then

error1(ea);

until ec=0;

readln;

repeat

winwwodab(ea);

str(a:4:2,f);

outtextxy(120,244,f);

gotoxy(16,18);

read(k);

val(k,b,ec);

if ec<>0 then

error1(ea);

until ec=0;

end;

procedure helpwin(ea:word);

{основа окна помощи}

begin

setfillstyle(1,3);

bar(350,100,590,380);

setcolor(0);

rectangle(353,103,587,377);

rectangle(355,105,585,375);

setcolor(14);

if (ea mod 2)=0 then

outtextxy(360,115,' ОКНО ПОМОЩИ')

else

outtextxy(360,115,' HELP WINDOW');

end;

procedure error1(ea:word);

begin

setcolor(15);

setfillstyle(1,12);

bar(140,210,490,280);

rectangle(145,215,485,275);

rectangle(147,217,483,273);

if (ea mod 2)=0 then

begin

outtextxy(150,227,' Ошибка! ');

outtextxy(150,237,' Вводимые параметр не число!! ');

outtextxy(150,250,' Проверьте значение и заново введите его.');

end

else

begin

outtextxy(150,227,' Error! ');

outtextxy(150,237,' The value you entered isn`t a quantity!!');

outtextxy(150,250,' Check it and put it in again. ');

end;

sound(600);

delay(4000);

nosound;

readln;

readln;

end;

procedure error(ea:word);

{Процедура ошибки}

begin

setcolor(15);

setfillstyle(1,12);

bar(140,210,490,260);

rectangle(145,215,485,255);

rectangle(147,217,483,253);

if (ea mod 2)=0 then

begin

outtextxy(150,227,' Ошибка!');

outtextxy(150,237,' Недостаток вводимых параметров!!');

end

else

begin

outtextxy(150,227,' Error!');

outtextxy(150,237,' Not all parameters are set!');

end;

sound(600);

delay(4000);

nosound;

readln;

end;

procedure newsctext(ea:word);

{Текст для процедуры newsc}

begin

if ea mod 2 =0 then

begin

settextstyle(0,0,1);

setcolor(15);

outtextxy(400,440,'Язык - Русский. ');

outtextxy(400,450,'Версия 1.0 Последнее издание');

outtextxy(400,460,'й Все права защищены.');

end

else

begin

settextstyle(0,0,1);

setcolor(15);

outtextxy(400,440,'Language - English.');

outtextxy(400,450,'Version 1.0 Final release.');

outtextxy(400,460,'й All rights reserved.');

end;

end;

procedure newsc(ea:word);

{Процедура обновления экрана}

begin

cleardevice;

setfillstyle(10,8);

floodfill(1,1,15);

setcolor(0);

setfillstyle(1,7);

bar(80,10,580,80);

rectangle(82,12,578,78);

rectangle(85,15,575,75);

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