Погрешность численного интегрирования

0.9992

99.9195

E-001 0.4137750613

E-001 0.2933317183

E-001 0.2032723690

E-001 0.4966040522

E-001 0.0896208493

E-001 0.4137750613

E-001 0.2933317183

E-001 0.2032723690

E-001 0.4966040522

E-001 0.0896208493

1.0000

E-015

1.0000

Else

1.0000

0.9996

E-004

1.0004

Интегрирование функций, заданных аналитически (формула прямоугольников, формула трапеций, формула Симпсона)

С геометрической точки зрения определенный интеграл

, (6.11)

есть площадь фигуры, ограниченная графиком функции и прямыми , (рис. 6.3).

Рис. 6.3

Разделим отрезок на N равных отрезков длиной , где

. (6.12)

Тогда координата правого конца i -го отрезка определяется по формуле

, (6.13)

где , .

Рис. 6.4. Метод левых прямоугольников	Рис. 6.5. Метод правых прямоугольников

Простейшая оценка площади под кривой может быть получена, как сумма площадей прямоугольников, одна из сторон которого совпадает с отрезком, а высота равна значению функции в точке (метод левых прямоугольников) (рис. 6.4) или в точке (метод правых прямоугольников) (рис. 6.5). (Погрешность вычисления значения интеграла на каждом шаге показана закрашенными фигурами.)

Значение определенного интеграла вычисляется по формулам

, (6.14)

(6.15)

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

Можно повысить точность вычисления определенного интеграла, если заменять реальную функцию на каждом интервале , , отрезком прямой, проходящей через точки с координатами ,(линейная интерполяция).

В этом случае фигура, ограниченная графиком функции и прямыми , , является трапецией. Искомый определенный интеграл определяется как сумма площадей всех трапеций:

. (6.16)

Более высокая точность вычисления интегралов обеспечивается при использовании параболической интерполяции (полиномом второй степени) по трем соседним точкам:

(6.17)

Более высокая точность вычисления интегралов обеспечивается при использовании параболической интерполяции (полиномом второй степени) по трем соседним точкам:

(6.17)

Для нахождения коэффициентов a, b, c полинома, проходящего через точки , , , нужно найти решение следующей системы линейных уравнений:

(6.18)

относительно неизвестных a, b, c.

Решив систему (6.18) относительно неизвестных a, b, c любым известным методом (например, Крамера), подставив найденные выражения в (6.17) и выполнив элементарные преобразования, получаем . (6.19)

Площадь под параболой на интервале находится посредством элементарного интегрирования (6.19):

, (6.20)

где . Искомый определенный интеграл находится как площадь всех параболических сегментов (формула Симпсона):

(6.21)

Обратите внимание на то обстоятельство, что в формуле Симпсона N должно быть четным числом.

Пример 6.2. Вычисление интеграла методом прямоугольников в пакете MATLAB:

>> f=inline('sin(x)'); % задание подынтегральной функции

>> Xmin=0;

>> Xmax=pi/2;

>> N=2001;

>> i=1:N;

>> dx=(Xmax-Xmin)/(N-1); % шаг интегрирования

>> x=Xmin:dx:Xmax; % вычисление координат узлов сетки

>> y=feval(f,x); % вычисление значений функции в узлах сетки

% вычисление интеграла по формуле правых прямоугольников

>> m=2:N;

>> y1(m-1)=y(m);

>> Fr=sum(y1)*dx

Fr =

>> Fr-1

ans =

% вычисление интеграла по формуле левых прямоугольников

>> m=1:N-1;

>> y1(m)=y(m);

>> Fl=sum(y1)*dx

Fl =

>> Fl-1

ans =

-3.9295e-004

% вычисление интеграла методом трапеций

>> s=0;

for i=2:N-1

s=s+y(i);

end;

Ft=(0.5*y(1)+s+0.5*y(N))*dx

Ft =

>> Ft-1

ans =

-5.1456e-008

% вычисление интеграла методом Симпсона

>> s=0;

for i=2:N-1

if i-2*ceil(i/2)==0

k=4;

k=2;

end;

s=s+k*y(i);

end;

Fs=(y(1)+s+y(N))*dx/3

Fs =

>> Fs-1

ans =

Для вычисления значений определенных интегралов в пакете MATLAB имеются следующие функции quad(), quadl(), trapz(), cumtrapz(), которые имеют следующий синтаксис.

>> q = quad(fun,a,b) % возвращает значение интеграла от функции

% fun на интервале [a,b], при вычислении

% используется адаптивный метод Симпсона.

>> q = quad(fun,a,b,tol) % возвращает значение интеграла от

% функции fun с заданной относительной

% погрешностью tol (по умолчанию tol=10^-3)

>> q = quad(fun,a,b,tol,trace) % возвращает значение интеграла от

% функции fun на интервале [a,b] на

% каждом шаге итерационного

% процесса

>> q = quad(fun,a,b,tol,trace,p1,p2,...) % возвращает значение

% интеграла от функции fun

% на на интервале [a,b] на

% каждом шаге

% итерационного процесса,

% p1, p2, … - параметры,

% передаваемые в функцию

% fun

>> [q,fcnt] = quadl(fun,a,b,...) % возвращает в переменную fcnt

% дополнительно к значению

% интеграла число выполненных

% итераций

Функция quadl() возвращает значения интеграла, используя для вычислений метод Лоббато (Lobbato). Синтаксис функции аналогичен синтаксису функции quad().

q = quadl(fun,a,b)

q = quadl(fun,a,b,tol)

q = quadl(fun,a,b,tol, 'trace ')

q = quadl(fun,a,b,tol, 'trace ',p1,p2,...)

[q,fcnt] = quadl(fun,a,b,...)

Пример 6.3. Вычисление интеграла с использованием функций quad.

>> q=quad('sin',0,pi/2,10^-4)

q =

>> q-1

ans =

-3.7216e-008

>> q=quad('sin',0,pi/2,10^-6,'trace');

>> q-1

ans =

-2.1269e-009

>> [q,fnct]=quad('sin',0,pi/2,10^-6,'trace');

>> 17

Функция trapz() вычисляет интеграл, используя метод трапеций. Синтаксис функции trapz():

Z = trapz(Y) % возвращает значение определенного интеграла, в

% предположении, что X=1:length(Y)

Z = trapz(X,Y) % возвращает значение интеграла

% на интервале [X(1),X(N)]

Z = trapz(...,dim) % интегрирует вектор Y, формируемый из чисел,

% расположенных в размерности dim

% многомерного массива

Функция cumtrapz() вычисляет интеграл, как функцию с переменным верхним пределом. Синтаксис функции cumtrapz() аналогичен синтаксису функции trapz().

Z = cumtrapz(Y)

Z = cumtrapz(X,Y)

Z = cumtrapz(... dim)

Пример 6.4. Вычисление определенного интеграла с использованием встроенной функции пакета MATLAB.

>> x=0:0.01:pi/2; % задание координат узловых точек

>> y=sin(x); % вычисление значений подынтегральной функции в

% узловых точках

>> trapz(y) % вычисление значения интеграла, в предположении о

% том, что шаг интегрирования равен единице

ans =

>> trapz(x,y) % вычисление значения интеграла на отрезке с

% шагом интегрирования 0.01

ans =

Пример 6.5. Вычисление интеграла с переменным верхним пределом

>> x=0:0.01:3*pi/2; % задание координат узловых точек

>> y=sin(x); % вычисление значений подынтегральной функции в

% узловых точках

>> z=cumtrapz(x,y); % вычисление значений интеграла с

% переменным верхним пределом в узловых

% точках

>> plot(x,y,x,z) % построение графиков подынтегральной функции и

% интеграла с переменным верхним пределом

% (рис. 6.6)

Рис. 6.6

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

. (6.22)

Тогда интеграл от данной функции на отрезке будет равен

. (6.23)

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

. (6.24)

Из (6.24) видно, что основной член погрешности на каждом отрезке имеет порядок или, в символической записи, . Поскольку полное число отрезков равно N, а , то полная погрешность метода левых прямоугольников по порядку величины равна . Аналогично можно показать, что погрешность метода правых прямоугольников также пропорциональна .

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

(6.25)

Заменив в (6.25) первый член выражением (6.23), значение функции в точке - разложением в ряд Тэйлора:

,

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

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

Полезно получить оценку погрешности вычисления интеграла от функции, зависящей от двух переменных, который с геометрической точки зрения представляет собой объем фигуры под поверхностью, заданной функцией . В прямоугольном приближении данный интеграл равен сумме объемов параллелепипедов с площадью основания и высотой, равной значению функции в одном из углов. Для определения погрешности разложим функцию в ряд Тейлора:

, (6.26)

где , - частные производные по соответствующим переменным.

Погрешность вычисления интеграла равна

. (6.27)

Подставив (6.26) в (6.27), выполнив интегрирование и приведя подобные, получаем, что член пропорциональный сокращается, а интеграл от дает . Интеграл от данного выражения по дает еще один множитель . Аналогичный вклад дает интеграл от члена, пропорционального . Так как погрешность порядок погрешности также составляет , то погрешность интегрирования по прямоугольнику , равна

. (6.28)

Из (6.28) видно, что погрешность интегрирования по одному параллелепипеду составляет . Так как имеется N параллелепипедов, полная погрешность по порядку величины равна . Однако в двумерном случае , поэтому полная погрешность . Напомним, что в одномерном случае полная погрешность метода прямоугольников .

Аналогичные оценки для двумерных обобщений формул трапеций и Симпсона показывают, что они соответственно равны и . Вообще можно показать, что если для одномерного случая погрешность составляет , то в d -мерном случае она равна .

Дата добавления: 2014-01-06