§ 2. Численное интегрирование по простейшим формулам

(((2*q-9)*q+11)*q-3)*r4[i]/12;

2: new1 := r2[i]+(q-1)*r3[i]+

((6*q-8)*q+11)*r4[i]/12;

end;

end.

{***Вторая интерполяционная формула Ньютона ***}

function new2 (key:integer;x,y,r1,r2,r3,r4:mas;

x1:real):real;

label 30;

var i,j,k : integer; q : real;

begin

i := 1;

j := 14;

if x1 > x[j] then goto 30;

{ **** Метод двоичного поиска интервала ****}

repeat

k := (i+j) div 2;

if x1< x[k] then j := k;

if x1>=x[k] then i := k;

until j <= i+1;

30:

q := (x1 - x[i]) / (x[2]-x[1]);

{ **** Выбор типа работы процедуры ****}

case key of

0: new2 := y[j]+q*(r1[j-1]+

(q+1)*(r2[j-2]/2.0+(q+2)*r3[j-3]/6.0));

1: new2 := r1[j-1]+0.5*(2*q+

1)*r2[j-2]+((3*q+6)*q+2)*r3[j-3]/6+

(((2*q+9)*q+11)*q+3)*r4[j-4]/12;

2: new2 := r2[j-2]+(q+

1)*r3[j-3]+((6*q+18)*q+11)*r4[j-4]/12;

end;

end.

Для проверки и тестирования процедур вы­чис­ля­лись зна­чения первой и второй про­из­вод­ных функ­ции при за­дан­ных значениях ар­гу­мен­та с ис­поль­зо­ва­­­нием пер­вой или вто­рой ин­тер­по­ля­­ци­он­ных фор­мул Нью­то­на. Расчеты конт­ро­­ли­ро­вались при со­с­тав­­ле­нии таблицы раз­ностей. Вы­чис­ле­ния вы­по­л­ня­лись с точ­ностью 10-5. Дан­­ные взяты из табл. 2.2. Ре­­зуль­та­ты работы про­це­дур при­водятся в табл. 3.1.

Таблица 3.1

Первая инт.формула Ньютона

Вторая инт.формула Ньютона

xi

y’

y’’

xi

y’

y’’

1.90800

-0.01067

-0.00043

2.24800

-0.01679

-0.00041

2.13500

-0.01449

-0.00042

2.35900

-0.01999

0.00025

§ 2. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ПО ПРОСТЕЙШИМ ФОРМУЛАМ

Пусть теперь некоторый конечный интервал [а, b] на оси Ох разбит на n подынтервалов [хi, хi+1], ко­то­рые в дальнейшем будем называть эле­мен­тар­ны­ми от­рез­ка­ми. Ясно, что без огра­ни­че­ния об­щ­нос­­ти можно по­ло­жить х0 = а; хn= b и х0 < <х1 < ... <хn. Через hi обозначим дли­­ну эле­мен­тар­но­го от­рез­ка (хi+1 - хi). Если заданный от­­ре­зок [а, b] разбит рав­но­мерно, то тогда hi будет по­сто­­ян­на для любой [а, b]. Пусть те­перь на [а, b] опре­де­ле­­на не­­которая функ­ция f(х). Пред­по­ло­жим, что не­об­хо­ди­мо най­ти приближение к определенному ин­­те­г­ра­­лу, которое обозначим . Оче­вид­но так­же, что ес­ли f(х) не­­прерывна [а, b], то тогда мож­­но пре­д­­­­­ста­вить как , где - ин­те­грал функ­­­ции f(х) на эле­ментарном от­рез­ке [хi+1, хi], т.е.:

.

Bсякая прос­­тая формула, ап­прок­си­ми­рующая от­­­дель­ный ин­теграл, на­зы­­вaется квад­ра­тур­ной. Составная квад­ра­­тур­ная формула - это фор­му­­­ла, да­ющая при­бли­же­ние ин­­­те­­гра­ла в ви­де суммы при­­ближений ин­те­­гра­ла­ми по дан­­ной квад­ра­тур­ной формуле.

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

Известны три раз­но­вид­­ности формул пря­мо­у­голь­­ни­ков: это формулы левых, правых и средних пря­­­мо­у­гольников. Все они основаны на ап­про­к­си­ма­­­ции каж­до­го ин­те­гра­ла площадью пря­мо­у­голь­­ника, од­ной из сто­рон которо­го яв­ляется hi , а вто­­рой - ли­бо значение функ­ции на левом кон­­це отрез­ка (рис.3.1, а), либо зна­че­ние функ­ции на пра­вом кон­це от­резка (рис. 3.1, б), либо зна­чение функ­­ции в сред­ней точке от­резка (рис. 3.1, в).

Квадратурные формулы, ап­прок­си­ми­ру­ю­щие , бу­дут иметь вид:

левых прямоугольников: = hi f (хi);

правых прямоугольников: = hi f (хi+1);

средних прямоугольников: = hi f (хi+1/2).

С учетом представления на элементарном от­­рез­ке со­став­ные квадpатурные формулы пря­мо­­у­голь­ни­ков мо­гут быть записаны так:

левых прямоугольников

;

правых прямоугольников

;

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

.

В формуле тра­пе­ций ис­­пользуются зна­­­че­ния функ­­­ции в кон­це­вых точ­ках эле­мен­тар­ных от­рез­ков. В этом случае ап­прок­си­­ми­­руется пло­ща­дь­ю тра­­­­пе­­­ции с ос­но­вани­я­­ми f(хi) и f(хi+1) и вы­­со­той Dx (рис. 3.2).

Тогда пло­щадь фигу­ры мо­­­жет быть опре­де­ле­на из фор­­­му­лы пло­ща­ди пря­­мо­у­голь­ной тра­пе­ции

Si = (fi + fi+1) hi /2.

Если теперь про­сум­ми­руем по­след­нюю фор­­му­лу по всем эле­мен­тар­ным от­рез­кам, то по­лу­чим с уче­том вы­пол­ненных эле­мен­тарных пре­об­­ра­зо­ваний сле­ду­ю­щее выражение:

Заметим, что при бесконечном умень­шении длин эле­­мен­тарных от­рез­ков формулы обоих типов (пря­мо­­у­голь­ни­ков и трапеций) сходятся к точ­но­му зна­чению ин­­те­гра­ла. Однако не яс­но, как быстро они сходятся ? По­пы­­та­ем­­­ся выяснить дан­ный во­прос, вос­­поль­зо­вав­шись раз­ло­же­ни­­ем фун­к­­­­ции в ряд Тейлора от­но­си­тель­но цен­тра эле­мен­тар­­ного от­резка [хi, хi+1]

f(x) = f(yi) + (x - yi) f (yi) + (x - yi)2 f ’’(yi)/2 +

+ (x - yi)3 f’’’(yi)/6 + (x-yi)4 f (IV)(yi)/24 + ...

затем, проинтегрировав по­лу­чен­ный ряд по каж­до­му из от­резков в пред­­положении, что ос­тав­ши­е­ся чле­ны ря­да на­мно­го меньше вы­­­пи­сан­ных, с учетом значений ко­эф­фи­ци­ен­тов ряда Тейлора на элементарном отрезке

получим

Член показывает ошибку фор­му­лы пря­­мо­угольников без учета членов более вы­со­кого по­ряд­ка. Ес­ли теперь подставим в фор­му­лу тра­пе­ций зна­чения функ­­ции в точках х = хi и х = хi+1, то по­лучим

Можно видеть, что ошибки формул сред­них пря­мо­у­голь­ни­­ков и трапеций одного по­ряд­ка, т.е. да­ют поч­ти оди­на­ко­­вую точность.

Примеров вы­чис­ле­ния определенных ин­те­г­­ра­лов по про­­с­тей­шим формулам (пря­мо­угольников и тра­пе­­ций) мож­­но при­вес­ти множество. Ни од­на БСП не об­­хо­дит­ся без этих про­цедур. Од­ним из примеров та­кой вы­чис­ли­тель­ной про­це­ду­­ры может слу­жить сле­­дующая:

procedure intlps (var s:real; k:integer);

begin

s := 0.0;

x := a;

for i := 0 to n do

begin

if k <>3 then x := a+ h*i

else x:=a+h/2+h*i;

f := func (a,b,x);

if k=4 then f := f/2.0;

if (i=0) and (k<>2) then s := s + f;

if (i=n) and ((k=2) or (k=4)) then s := s + f;

if (i<>0) and (i<>n) then

if k<>4 then s := s+f

else s := s+f*2;

end; s := s * h;

end.

Эта процедура обращается к процедуре-функ­ции func для вычисления подынтегральной функ­ции в заданных

71

Соседние файлы в папке GLAVA3_1