Погрешность численного интегрирования
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-мерном случае она равна 
.