Версия от 18.10.10
Методы численного интегрирования
Первообразной функцией для функции f(x) называется функция
F(x) такая, что F'(x)=f(x).
Если функция f(x) непрерывна на отрезке [a, b] и известна её первообразная функция F(x), то значение определённого интеграла от этой функции в пределах от a до b может быть вычислено по формуле Ньютона-Лейбница
b |
|
∫ f ( x )dx = F( b ) − F( a ) |
(1) |
b
Во многих случаях первообразная функция F(x) не может быть найдена с помощью известных средств или является слишком сложной, делая вычисление интеграла по формуле (1) чрезвычайно затруднительным или даже практически невыполнимым.
Если же подынтегральная функция f(x) задана таблично, то понятие первообразной теряет смысл.
Численное вычисление однократного интеграла называется
механической квадратурой.
Обычный приём механической квадратуры состоит в том, что функцию f(x) на рассматриваемом интервале [a, b] заменяют интерполирующей или аппроксимирующей функцией φ(x) простого вида (например, полиномом), а затем полагают
∫b |
f ( x )dx ∫b ϕ( x )dx |
(2) |
a |
a |
|
|
|
b |
Функция φ(x) должна быть такой, что |
∫ϕ( x )dx вычисляется |
|
непосредственно [1]. |
a |
|
Квадратурные формулы Ньютона-Котеса
Формулы Ньютона-Котеса требуют чтобы значения аргумента х
винтервале интегрирования [a, b] были заданы с постоянным шагом.
Втеоретическом курсе высшей математики интеграл определяется как предел
b |
f ( x )dx = |
lim |
n |
f ( x ) |
x + R = |
n |
S |
|
+ R |
|
∫ |
∑ |
∑ |
|
(3) |
||||||
|
max xi →0 |
i |
i |
|
i |
|
|
|||
a |
|
|
i=0 |
|
|
i=0 |
|
|
|
|
где R - погрешность вычисления интеграла.
Простейшие квадратурные формулы можно получить исходя из простых наглядных соображений.
1 |
Любимов Е.Б. |
|
Версия от 18.10.10
Разобьём интервал интегрирования [a, b] на n равных отрезков
(рис. 1)
xi=a+i*h, (i = 0, 1, 2, . . ., n).
Вэтом случае значение xn=b, а значение h=(b-a)/n.
Впределах интервала [xi , xi+1] можно считать f(х)≈f(ξ)≈const, а точка ξ Є [xi , xi+1]. Этот подход определяет множество методов, определяемых как методы прямоугольников. Наиболее часто используются методы левых, правых и средних прямоугольников, в которых положение точки ξ определяется соответственно левой, правой границей интервала [xi , xi+1] или в точке, лежащей в середине
отрезка (xi +xi+1)/2.
Рис.1. Графическое представление интегрирования методом левых прямоугольников
Рассмотрим простой пример вычисления значения интеграла, в котором используются вышеперечисленные методы численного интегрирования.
Пусть требуется вычислить значение интеграла
∫1 |
dx |
|
1 + x |
||
0 |
Положим значение n=6. В этом случае значение h=(1-0)/6=0,166(6).
2 |
Любимов Е.Б. |
|
Версия от 18.10.10
Вычислим таблицу значений подынтегральной функции в узлах xi=x0+i*h и в точках, определяющих середины подынтервалов
(xi +xi+1)/2.
h= |
0,16667 |
|
S лев |
|
|
|
S cp |
i |
x i |
f(x i ) |
0 |
S прав |
(x i +x i+1 )/2 |
f(x cp ) |
0 |
0 |
0 |
1 |
0,16667 |
0 |
0,0833 |
0,92308 |
0,15385 |
1 |
0,16667 |
0,85714 |
0,30952 |
0,14286 |
0,25 |
0,8 |
0,28718 |
2 |
0,33333 |
0,75 |
0,43452 |
0,26786 |
0,4167 |
0,70588 |
0,40483 |
3 |
0,5 |
0,66667 |
0,54563 |
0,37897 |
0,5833 |
0,63158 |
0,51009 |
4 |
0,66667 |
0,6 |
0,64563 |
0,47897 |
0,75 |
0,57143 |
0,60533 |
5 |
0,83333 |
0,54545 |
0,73654 |
0,56988 |
0,9167 |
0,52174 |
0,69228 |
6 |
1 |
0,5 |
|
0,65321 |
|
|
0,69228 |
S лев = |
0,73654 |
|
S прав = |
0,65321 |
|
S cp = |
0,69228 |
I=ln(2)= |
0,693147 |
- точное значение интеграла |
|
|
|
||
Для оценки точности вычисленного значения интеграла осуществляются повторные вычисления интегралов при значениях n, увеличенных в два раза.
h= |
0,08333 |
|
S лев |
|
|
|
S cp |
i |
x i |
f(x i ) |
0 |
S прав |
(x i +x i+1 )/2 |
f(x cp ) |
0 |
0 |
0 |
1 |
0,08333 |
0 |
0,0417 |
0,9600 |
0,0800 |
1 |
0,08333 |
0,92308 |
0,16026 |
0,076923 |
0,1250 |
0,8889 |
0,1541 |
2 |
0,16667 |
0,85714 |
0,23168 |
0,148352 |
0,2083 |
0,8276 |
0,2230 |
3 |
0,25000 |
0,8 |
0,29835 |
0,215018 |
0,2917 |
0,7742 |
0,2876 |
4 |
0,33333 |
0,75 |
0,36085 |
0,277518 |
0,3750 |
0,7273 |
0,3482 |
5 |
0,41667 |
0,70588 |
0,41968 |
0,336342 |
0,4583 |
0,6857 |
0,4053 |
6 |
0,50000 |
0,66667 |
0,47523 |
0,391897 |
0,5417 |
0,6486 |
0,4594 |
7 |
0,58333 |
0,63158 |
0,52786 |
0,444529 |
0,6250 |
0,6154 |
0,5106 |
8 |
0,66667 |
0,6 |
0,57786 |
0,494529 |
0,7083 |
0,5854 |
0,5594 |
9 |
0,75000 |
0,57143 |
0,62548 |
0,542148 |
0,7917 |
0,5581 |
0,6059 |
10 |
0,83333 |
0,54545 |
0,67094 |
0,587603 |
0,8750 |
0,5333 |
0,6504 |
11 |
0,91667 |
0,52174 |
0,71441 |
0,631081 |
0,9583 |
0,5106 |
0,6929 |
12 |
1,00000 |
0,5 |
|
0,672747 |
|
|
|
|
|
|
|
|
|
|
|
S лев = |
0,714414 |
|
S прав = |
0,672747 |
|
S cp = |
0,692930 |
Метод прямоугольников может быть реализован в трёх вариантах: левых прямоугольников, правых прямоугольников и средних прямоугольников. В формулах метода прямоугольников одна
3 |
Любимов Е.Б. |
|
Версия от 18.10.10
сторона каждого прямоугольника h определяется по формуле (3), а значение второй стороны в зависимости от метода равно:
f(xi) для i=0, 1, 2,...,n-1 в формуле левых прямоугольников
|
b |
n−1 |
|
|
|
n−1 |
|
|
|
∫ f ( x )dx = ∑f ( xi )h + RЛП = h∑f ( a +i • h ) + RЛП |
(4) |
||||||
|
a |
i=0 |
|
|
|
i=0 |
|
|
∫b |
f(xi) для i= 1, 2,...,n в формуле правых прямоугольников |
|
||||||
f ( x )dx = ∑n |
f ( xi )h + RПрП = h∑n |
f ( a +i • h ) + RПрП |
(5) |
|||||
a |
i=1 |
|
|
i=1 |
|
|
|
|
|
и f([xi + xi+1] / 2) для i=0, 1, 2,...,n-1 в формуле средних |
|
||||||
|
прямоугольников |
|
|
|
|
|
|
|
|
|
b |
|
|
xi |
+ xi+1 |
|
|
|
|
∫ f ( x )dx = h∑n |
f ( |
)h + RСрП |
(6) |
|||
|
|
|
2 |
|||||
|
|
a |
i=1 |
|
|
|
|
|
Обратите внимание на индексы в определении значения аргумента xi в обращениях к функциям, используемых для вычисления площадей элементарных прямоугольников и на то, что во всех вариантах метода суммируются площади (n-1) - го прямоугольника.
Как видно из приведённого выше примера, наиболее точное значение интеграла получено при использовании метода средних прямоугольников.
Метод трапеций
Рис. 2. Элементарная трапеция, площадь которой вычисляется в методе трапеций
4 |
Любимов Е.Б. |
|
Версия от 18.10.10
Площадь элементарной фигуры моделируется площадью трапеции, у которой высота равна h, а параллельные стороны равны соответственно f(xi) и f(xi+1). Площадь этой элементарной трапеции вычисляется по формуле
S |
i |
= h |
f (xi ) + f (xi+1 ) |
+ R , (i = 0,1,2,..., n −1) |
(7) |
|
|||||
|
|
2 |
i |
|
|
|
|
|
|
|
где Ri площадь, сегмента залитого на рис.2 серым цветом.
Значение интеграла, вычисляемое по формуле метода трапеций
равно |
b |
n−1 |
n−1 |
|
|
|
|||
|
∫f (x)dx =∑Si +∑Ri =S +R |
(8) |
||
|
a |
i=0 |
i=0 |
|
Выполнив дополнительные алгебраические преобразования, получим окончательный вид формулы метода трапеций:
b |
n−1 |
h |
|
n−1 |
|
|
∫f ( x )dx ∑Si = |
[f ( a )+ f (b)]+h∑f ( xi ), |
(9) |
||||
|
||||||
a |
i=0 |
2 |
i=1 |
|
||
где xi=a+ih ( i=1,2,…,n-1).
h= |
0,16667 |
|
Sтрап |
i |
x i |
f(x i ) 0,00000 |
|
00 1 0,15476
10,16667 0,85714 0,28869
20,33333 0,75 0,40675
30,5 0,66667 0,51230
4 |
0,66667 |
0,6 |
0,60776 |
5 |
0,83333 |
0,54545 |
0,69488 |
6 |
1 |
0,5 |
|
|
|
|
Sтрап |
h= |
0,08333 |
|
|
i |
x i |
f(x i ) |
0,00000 |
0 |
0,00000 |
1 |
0,08013 |
1 |
0,08333 |
0,92308 |
0,15430 |
2 |
0,16667 |
0,85714 |
0,22335 |
3 |
0,25000 |
0,8 |
0,28793 |
4 |
0,33333 |
0,75 |
0,34860 |
5 |
0,41667 |
0,70588 |
0,40579 |
6 |
0,50000 |
0,66667 |
0,45988 |
7 |
0,58333 |
0,63158 |
0,51120 |
8 |
0,66667 |
0,6 |
0,56001 |
9 |
0,75000 |
0,57143 |
0,60654 |
10 |
0,83333 |
0,54545 |
0,65101 |
11 |
0,91667 |
0,52174 |
0,69358 |
12 |
1,00000 |
0,5 |
|
5 |
Любимов Е.Б. |
|