Приближенное вычисление определенного интеграла при помощи квадратурной формулы Чебышева

МИНИСТЕРСТВО ОБРАЗОВАНИЯ УКРАИНЫ

ГОСУДАРСТВЕННЫЙ ХИМИКОТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ

КАФЕДРА ВЫЧИСЛИТЕЛЬНОЙ ТЕХНИКИ

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

Приближенное вычисление определенного интеграла при помощи квадратурной формулы Чебышева

Выполнил:

студент 2-го курса

Поляков Е.В.

Днепропетровск,

2000

1. Общая постановка и анализ задачи

1.1 Введение

Требуется найти определенный интеграл:

по квадратурной формуле Чебышева.

Рассмотрим, что представляет собой вообще квадратурная формула, и как можно с ее помощью вычислить приближенно интеграл.

Известно, что определенный интеграл функции типа численно представляет собой площадь криволинейной трапеции ограниченной кривыми x=0, y=a, y=b и y= (Рис.1).

Рис. 1. - Криволинейная трапеция

Если f(x) непрерывна на отрезке [a, b], и известна ее первообразная F(x), то определенный интеграл от этой функции в пределах от а до b может быть вычислен по, известной всем, формуле Ньютона-Лейбница:

- F(a),

где F'(x) = f(x).

Однако во многих случаях F(x) не может быть найдена, или первообразная получается очень сложной для вычисления.

Кроме того, функция часто задается таблично. Поэтому большое значение приобретает приближенное и в первую очередь численное интегрирование.

Задача численного интегрирования состоит в нахождении приближенного значения интеграла по заданным или вычисленным значениям подинтегральной функции f(x) в некоторых точках (узлах) отрезка [a, b].

Численное определение однократного интеграла называется механической квадратурой, а соответствующие формулы численного интегрирования - квадратурными.

Заменяя подинтегральную функцию каким-либо интерполционным многочленом, мы получим квадратурные формулы вида:

где xk - выбранные узлы интерполяции; Ak - коэффициенты, зависящие только от выбора узлов, но не от вида функции (k=0,1,2,........, n); R - остаточный член, или погрешность квадратурной формулы.

Отбрасывая остаточный член R, мы совершаем погрешность усечения. При расчете к ней добавляются еще различные погрешности округления.

Разобьем отрезок интегрирования [a, b] на n равных частей системой точек:

xi= xo+ i..h; ( i = 0,1,2,......,n)

xo=a; xn=b;

h=(b-a)/n

и вычислим подинтегральную функцию в полученных узлах:

yi= f(xi) ; ( i = 0,1,2,......,n).

1.2 Вывод формул численного интегрирования с использованием интерполяционного полинома Лагранжа

Пусть для y=f(x) известны в n+1 точках X0,X1,X2..Xn промежутка [a,b] соответствующие значения f(xi)=yi (i=0,1,2..n). Требуется приближенно найти:

По заданным значениям Yi построим полином Лагранжа. Заменим f(x) полиномом Ln(x). Тогда

,

где Rn(f) - ошибка квадратурной формулы.

Отсюда, воспользовавшись выражением для Ln(x), получаем приближенную квадратурную формулу:

Для вычисления коэффициентов Аi заметим что:

1. коэффициенты Ai при данном расположении узлов не зависит от выбора функции f(x);

2. для полинома степени n последняя формула точная.

Полаая y=xk (k=0,1,2..,n), получим линейную систему из n+1 уравнений:

,

где

, (k=0,1,..,n),

из которой можно определить коэффициенты А0,А1,..,АN.

Определитель системы есть определитель Вандермонда:

Заметим, что при применении этого метода фактическое построение полинома Лагранжа Ln(x) является излишним. Простой метод подсчета погрешности квадратурных формул разработан С.М. Никольским.

Теперь рассмотрим несколько простейших квадратурных формул:

1.3 Формула трапеций и средних прямоугольников

Заменим дугу АВ стягивающей ее хордой, получим прямолинейную трапецию аАВb, площадь которой примем за приближенное значение интеграла.

Рис 2. - Криволинейная трапеция

Рис. 3. - Метод трапеций



Рис. 4. - Метод средних прямоугольников

По методам трапеций и средних прямоугольников (рис. 3, рис. 4) соответственно интеграл равен сумме площадей прямоугольных трапеций, где основание трапеции какая-либо малая величина (точность), и сумма площадей прямоугольников, где основание прямоугольника какая-либо малая величина (точность), а высота определяется по точке пересечения верхнего основания прямоугольника, которое график функции должен пересекать в середине. Соответственно получаем формулы площадей:

- для метода трапеций:

,

- для метода средних прямоугольников:

.

1.4 Общая формула Симпсона (параболическая формула)

Пусть n=2m есть четное число и yi=f(xi) (i=0,1,2...n) - значения функции y=f(x) для равноотстоящих точек а=x0,x1, ... ,xn=b с шагом

.

Применив формулу Симпсона к каждому удвоенному промежутку [x0,x2], [x2,x4] ... [x2m-2,x2m] длины 2h и введя обозначения

1=y1+y2+ ... +y2m-1;

2=y2+y4+ ... +y2m,

получим обобщенную формулу Симпсона:

.

Остаточный член формулы Симпсона в общем виде:

,

где k=(x2к-2,x2к).

1.5 Квадратурная формула Чебышева

Рассмотрим квадратурную формулу вида:

Функцию f(x) будем искать в виде, когда f(x) многочлен: f(x)=ao+a1x+...+anxn . Проинтегрировав, преобразовав и подставив значения многочлена в узлах:

f(x1)=a0+a1x1+a2x12+a3x13+...+anx1n;

f(x2)=a0+a1x2+a2x22+a3x23+...+anx2n;

f(x3)=a0+a1x3+a2x32+a3x33+...+anx3n;

. . . . . . . . . . . . . . . .

f(xn)=a0+a1xn+a2xn2+a3xn3+...+anxnn.

Получим формулу Чебышева:

Значения х1, х2,..,хn для различных n приведены в таблице.

Значения х1,х2,..,хn для различных n

n	I	ti	n	i	ti
2	1;2	0,577350	6	1;6	0,866247
3	1;3	0,707107		2;5	0,422519
	2	0		3;4	0,266635
4	1;4	0,794654	7	1;7	0,883862
	2;3	0,187592		2;6	0,529657
5	1;5	0,832498		3;5	0,321912
	2;4	0,374541		4	0
	3	0

2. Решение контрольного примера

Здесь a=0 ; b= ; при n=5.

f(x)=sin(x);

;

.

i	xi	yi
1	0,131489	0,131118
2	0,490985	0,471494
3	0,785	0,706825
4	0,509015	0,487317
5	0,868511	0,763367

x1=/4+/4*t1=/4+/4(-0,832498)=0,131489;

x2=/4+/4*t2=/4+/4(-0,374341)=0,490985;

x3=/4+/4*t3=/4+/4*0=0,785;

x4=1-x2=1-0,490985 = 0,509015;

x5=1-x1=1-0,131489=0,868511;

y1=sin(x1)=sin(0,131489)=0,131118;

y2=sin(x2)=sin(0,490985)=0,471494;

y3=sin(x3)=sin(0,785)=0,706825;

y4=sin(x4)=sin(0,509015)=0,487317;

y5=sin(x5)=sin(0,868511)=0,763367;

;

I=/10(0,131118+0,471494+0,706825+0,487317+0,763367)=

=/10*2,560121=0,8038779.

3. Описание программы Integral. pas. Алгоритм

Процедура VVOD - заполняет массив, содержащий в себе аргументы xi.

Процедура FORM - используя массив, содержащий аргументы xi заполняет массив yi

Процедура CHEB - используя массивы xi и yi, высчитывает по квадратурной формуле Чебышева приближенное значение интеграла.

Процедура TABL - это подпрограмма, осуществляющая вывод таблицы узлов (аргумент - функция).

При запуске программы нужно ввести границы интегрирования.

После ввода границ интегрирования используется процедура VVOD, а затем высчитывается, и выводиться на экран шаг табулирования функции h.

После этого используем процедуры FORM и CHEB .

Получив результат, выводим таблицу (процедура TABL ) и интеграл.

4. Заключение и выводы

Таким образом, очевидно, что при вычислении определенных интегралов с помощью квадратурных формул, а в частности по формуле Чебышева не дает нам точного значения, а только приближенное.

Чтобы максимально приблизиться к достоверному значению интеграла нужно уметь правильно выбрать метод и формулу, по которой будет вестись расчет. Так же очень важно то, какой будет взят шаг интегрирования.

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

Листинг программы

Программа написана на языке Tubro Pascal 7.0 для MS-DOS. Ниже приведен ее листинг:

program integral;

uses crt;

const n=5;

k=-0.832498;

l=-0.374541;

z=0.0;

type aa=array[1..n] of real;

var x,y:aa;

a,b,h,ich:real;

{заполнение х-сов в массив х[5]}

procedure vvod(var a,b:real;var c:aa);

var i:integer;

t:aa;

Begin

t[1]:=k;

t[2]:=l;

t[3]:=z;

t[4]:=l;

t[5]:=k;

for i:=1 to n-1 do

c[i]:=((b+a)/2+(b-a)/2*t[i]);

for i:=n-1 to n do

c[i]:=1 - c[n+1-i];

end;

{заполнение y-ков в массиве у[5]}

procedure form(var x:aa; var y:aa);

var i:integer;

Begin

for i:=1 to n do

y[i]:=sin(x[i]); {функция}

end;

{процедура для расчета интеграла по квадратурной формуле Чебышева}

procedure cheb(var y:aa;var ich:real);

var i:integer;

Begin

ich:=0;

for i:=1 to n do

ich:=ich+y[i]*h;

end;

{процедура вывода таблицы}

procedure tabl;

var i:integer;

Begin

writeln(' ___________________________________ ');

writeln('| i | t| x|y |');

writeln(' ___________________________________ ');

writeln('| 1 |',k:9:6,'|',x[1]:9:6,' |',y[1]:9:6,'|');

writeln('| 2 |',l:9:6,'|',x[2]:9:6,' |',y[2]:9:6,'|');

writeln('| 3 |',z:9:6,'|',x[3]:9:6,' |',y[3]:9:6,'|');

writeln('| 4 |',l:9:6,'|',x[4]:9:6,' |',y[4]:9:6,'|');

writeln('| 5 |',k:9:6,'|',x[5]:9:6,' |',y[5]:9:6,'|');

writeln(' ___________________________________ ');

end;

Begin

clrscr;

writeln('ПРОГРАММА ДЛЯ ВЫЧИСЛЕНИЯ ОПРЕДЕЛЕННОГО ИНТЕГРАЛА');

writeln;

writeln('Введите границы интегрирования a,b:');

readln(a,b);

vvod(a,b,x);

h:=(b-a)/n;

writeln('h=',h:9:6);

form(x,y);

cheb(y,ich);

tabl;

writeln('I=',ich:8:6);

end.

Вывод результата:

ПРОГРАММА ДЛЯ ВЫЧИСЛЕНИЯ ОПРЕДЕЛЕННОГО ИНТЕГРАЛА

Введите границы интегрирования a,b:

0 1.5708

h= 0.314160

____________________________

| i | t | x | y |

____________________________

| 1 |-0.832498| 0.131556 | 0.131177|

| 2 |-0.374541| 0.491235 | 0.471716|

| 3 | 0.000000| 0.785400 | 0.707108|

| 4 |-0.374541| 0.508765 | 0.487099|

| 5 |-0.832498| 0.868444 | 0.763325|

____________________________

I=0.804383

Список литературы

1. Ракитин Т.А., Первушин В.А. “Практическое руководство по численным методам с приложением программ на языке Basic”.

2. Крылов В.И. “Приближенные вычисления интегралов” - М.: Физмат.

3. Демидович и Марон “Основы вычислительной математики”.

4. Копченова и Марон “Вычислительная математика в примерах и задачах”.

5. Вольвачев А.Н., Крисевич В.С. Программирование на языке Паскаль для ПЭВМ ЕС. Минск: 1989 г.

6. Зуев Е.А. Язык программирования Turbo Pascal. М.1992 г.

7. Скляров В.А. Знакомьтесь: Паскаль. М. 1988 г.

8. Материалы сайта http://www.ed.vseved.ru/.