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

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

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

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

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

Лицей №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) задана в уздах интерполяции:

(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

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

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;

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.');

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