Численное интегрирование
к.ф.-м.н. Уткин Павел Сергеевич 1
e-mail: [email protected], [email protected](926) 2766560
Данная лекция доступна по адресуhttp://mipt.ru/education/chair/computational_mathematics/study/materials/compmath/lectures/
4 октября 2014, МФТИ, Долгопрудный
1Конспект Ивана Цыбулина, email: [email protected]
Численное интегрирование Введение
Постановка задачи
Для заданной функции f (x) вычислить значение определенногоинтеграла
I =
∫baf (x)dx .
Желательно, чтобы метод численного интегрирования обладалследующими свойствами:• Универсальность. Функция f (x) может быть задана в виде«черного ящика», как способ вычисления f (x) по данному x .
• Экономичность. Количество вычислений функции f (x) повозможности должно быть сведено к минимуму.
• Хорошая обусловленность. Неустранимые погрешности ∆f взначениях f (x) не должны приводить к значительной итоговойошибке ∆I .
Уткин П.С. Интегрирование 2 / 45
Численное интегрирование Введение
Приложения
Численное интегрирование может применяться для• интегрирования функций, известных только в некоторых точках,например, полученных в результате измерений;
• интегрирования сложных выражений, не имеющих элементарныхпервообразных, либо имеющих слишком громоздкие выражениядля них;
• построения методов численного решения уравнений вобыкновенных и частных производных (методы конечныхэлементов, интегро-интерполяционные методы).
Уткин П.С. Интегрирование 3 / 45
Численное интегрирование Простейшие методы
Методы, основанные на определении интеграла
Вспомним, как определяется определенный интеграл Римана.Рассмотрим на отрезке интегрирования [a, b] сеткуa = x0, x1, . . . , xN = b. Возьмем от каждого отрезка по однойточке-представителю ξi ∈ [xi−1, xi ] и составим интегральную сумму
SN =
N∑i=1
f (ξi )∆xi
Значение интеграла определяется как предел при max∆xi → 0значений интегральных сумм:∫b
af (x)dx = lim
max∆xi→0
N∑i=1
f (ξi )∆xi
Уткин П.С. Интегрирование 4 / 45
Численное интегрирование Простейшие методы
Формула прямоугольников
Возьмем в качестве приближенного значения интеграла значениенекоторой интегральной суммы:∫b
af (x)dx ≈ SN =
N∑i=1
f (ξi )∆xi .
Возьмем в качестве ξi середину отрезка [xi−1, xi ]:∫baf (x)dx ≈
N∑i=1
f
(xi−1 + xi
2
)∆xi .
Методы численного интегрирования называют квадратурнымиформулами или просто квадратурами. Данный метод называетсяформулой средней точки или формулой средних прямоугольников.
Уткин П.С. Интегрирование 5 / 45
Численное интегрирование Простейшие методы
Формула средних прямоугольников
∫baf (x)dx ≈
N∑i=1
f
(xi−1 + xi
2
)∆xi .
Уткин П.С. Интегрирование 6 / 45
Численное интегрирование Простейшие методы
Формулы односторонних прямоугольников
Выбор в качестве ξi средней точки интервала не принципиален, можновзять, например, левый или правый конец интервала.Соответствующие формулы называются формулами левых и правыхпрямоугольников.
∫baf (x)dx ≈
N∑i=1
f (xi−1)∆xi .
∫baf (x)dx ≈
N∑i=1
f (xi )∆xi .
Уткин П.С. Интегрирование 7 / 45
Численное интегрирование Простейшие методы
Составные квадратурные формулы
Построенные формулы прямоугольников являются составнымиквадратурными формулами. Квадратурная формула называетсясоставной, если является результатом применения некоторойэлементарной квадратурной формулы к каждому интервалу [xi−1, xi ].
Элементарная квадратура Составная квадратура∫baf (x)dx ≈ (b − a)f
(a + b
2
) ∫baf (x)dx ≈
N∑i=1
f
(xi−1 + xi
2
)∆xi∫b
af (x)dx ≈ (b − a)f (a)
∫baf (x)dx ≈
N∑i=1
f (xi−1)∆xi∫baf (x)dx ≈ (b − a)f (b)
∫baf (x)dx ≈
N∑i=1
f (xi )∆xi
Будем дальше строить только элементарные квадратурные формулы,составные получаются из них тривиально.
Уткин П.С. Интегрирование 8 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Приближение подынтегральной функции
Пусть подынтегральная функция f (x) хорошо приближается некоторойпросто интегрируемой функцией P(x). Тогда∫b
af (x)dx ≈
∫baP(x)dx .
В качестве функции P(x) могут выступать• многочлены;• тригонометрические многочлены;• экспоненциальные многочлены и тому подобные.
Остановимся на случае, когда P(x) — многочлен. Тогда задачаприближения f (x) с помощью функции P(x) может быть решена,например, с помощью алгебраической интерполяции.
Уткин П.С. Интегрирование 9 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Интерполяционные квадратурные формулы
Введем на отрезке [a, b] некоторую сетку
Ωs = a 6 x0 < x1 < · · · < xs 6 b .
Построим на этой сетке интерполяционный многочлен в формеЛагранжа:
PL(x) =s∑
k=0
ck(x)f (xk),
где ck(x) — это базисные интерполяционные многочлены Лагранжа:
ck(x) =s∏
i=0i 6=k
x − xixk − xi
.
Интегрируя PL(x) по отрезку [a, b], имеем:∫baPL(x)dx = (b − a)
s∑k=0
wk f (xk), wk ≡1
b − a
∫back(x)dx
Уткин П.С. Интегрирование 10 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Интерполяционные квадратурные формулы
∫baf (x)dx ≈
∫baPL(x)dx = (b−a)
s∑k=0
wk f (xk), wk ≡1
b − a
∫back(x)dx
Утверждение 5.1.Величины wk не зависят от конкретных значений a, b, xi , а зависятлишь от относительного расположения узлов xi .
Доказательство.
Сделаем замену переменных x = a+b2 + ξb−a
2 , ξ ∈ [−1, 1]:
wk =1
b − a
∫ba
s∏i=0i 6=k
x − xixk − xi
dx =12
∫1−1
s∏i=0i 6=k
ξ− ξiξk − ξi
dξ
Уткин П.С. Интегрирование 11 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Интерполяционные квадратурные формулы
Вид квадратурной формулы∫baf (x)dx ≈ (b − a)
s∑k=0
wk f (xk)
является универсальным, и мог быть написан из общих соображенийлинейности квадратурной формулы по значениям подынтегральнойфункции (по аналогии с линейностью самого интеграла).Величины wk называются весами квадратурной формулы. Веса независят от конкретного отрезка интегрирования, и могут бытьвычислены на некотором «стандартном» отрезке. Обычно используютотрезок [0, 1] или [−1, 1].Соответственно xk называются узлами квадратурной формулы. Онитакже обычно приводятся для стандартного отрезка, а дляконкретного отрезка [a, b] они получаются линейным преобразованием.
Уткин П.С. Интегрирование 12 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Формулы Ньютона-Котеса
Изучим случай, когда Ω — равномерная сетка:
xi = a +b − a
si , i = 0, . . . , s.
Интерполяционные квадратурные формулы на такой сетке называютсяформулами Ньютона-Котеса. Некоторые из них имеют свои названия.
Название Узлы Веса Вид
Трапеций a, b 1/2 (b − a)f (a) + f (b)
2
Симпсонаa, b 1/6
(b − a)f (a) + 4f (a+b
2 ) + f (b)
6a+b2 2/3
Правило 3/8a, b 1/8 f (a) + 3f (2a+b
3 ) + 3f (a+2b3 ) + f (b)
8/(b − a)2a+b3 , a+2b
3 3/8
Уткин П.С. Интегрирование 13 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Формула трапеций
Для одного отрезка:∫baf (x)dx ≈ (b − a)
f (a) + f (b)
2
Составная формула:∫baf (x)dx ≈
N∑i=1
f (xi−1) + f (xi )
2∆xi
Уткин П.С. Интегрирование 14 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Формула Симпсона
Для одного отрезка:∫baf (x)dx ≈ (b − a)
f (a) + 4f (a+b2 ) + f (b)
6
Составная формула:∫baf (x)dx ≈
N∑i=1
f (xi−1) + 4f ( xi−1+xi2 ) + f (xi )
6∆xi
Уткин П.С. Интегрирование 15 / 45
Численное интегрирование Интерполяционные квадратурные формулы
Формула «3/8»
Для одного отрезка: (b − a)f (a) + 3f (2a+b
3 ) + 3f (a+2b3 ) + f (b)
8
Составная формула:N∑i=1
f (xi−1) + 3f (2xi−1+xi3 ) + 3f ( xi−1+2xi
3 ) + f (xi )
8∆xi
Уткин П.С. Интегрирование 16 / 45
Численное интегрирование Погрешность квадратурных формул
Что точнее?
Попробуйте угадать, какая из двух квадратурных формул точнее:
O x
y
O x
y
Метод средней точки
Ошибка ∼ 1%
Метод трапеций
Ошибка ∼ 2%
Уткин П.С. Интегрирование 17 / 45
Численное интегрирование Погрешность квадратурных формул
Что точнее?
Попробуйте угадать, какая из двух квадратурных формул точнее:
O x
y
O x
y
Метод средней точкиОшибка ∼ 1%
Метод трапецийОшибка ∼ 2%
Уткин П.С. Интегрирование 17 / 45
Численное интегрирование Погрешность квадратурных формул
Степень квадратурной формулы
Будем говорить, что квадратурная формула имеет алгебраическуюстепень точности m, если она точно интегрирует все многочленыстепени не более m, а некоторые многочлены степени m + 1 уже нет.Например, формула средней точки имеет алгебраическую степень 1:она точна для всех линейных функций∫b
a(αx + β)dx = (b − a)β+ α
b2 − a2
2= (b − a)
(αa + b
2+ β
)и не точна для 3x2∫b
a3x2dx = b3 − a3 = (b − a)(a2 + ab + b2) 6= (b − a)
34(a2 + 2ab + b2)
Уткин П.С. Интегрирование 18 / 45
Численное интегрирование Погрешность квадратурных формул
Условия заданной степени
Утверждение 5.2.Для того, чтобы квадратурная формула имела алгебраическуюстепень m необходимо и достаточно, чтобы она точно интегрировалафункции 1, x , x2, . . . , xm.
Доказательство.
Достаточность. Пусть Pm(x) = α0 + α1x + · · ·+ αmxm
(b − a)s∑
k=0
wkPm(x) = (b − a)s∑
k=0
wk
m∑r=0
αrxr =
=
m∑r=0
αr (b − a)s∑
k=0
wkxr =
m∑r=0
αr
∫bax rdx =
∫baPm(x)dx
Необходимость очевидна.
Уткин П.С. Интегрирование 19 / 45
Численное интегрирование Погрешность квадратурных формул
Условия заданной степени
Условия из утверждения 5.2 позволяют строить квадратурныеформулы степени m, решая систему уравнений (для простоты взятотрезок [−1, 1]):
2s∑
k=0
wk =
∫1−1
dx = 2
2s∑
k=0
xkwk =
∫1−1
xdx = 0
2s∑
k=0
x2kwk =
∫1−1
x2dx =23
...
2s∑
k=0
xmk wk =
∫1−1
xmdx =1+ (−1)m
m + 1
Уткин П.С. Интегрирование 20 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член квадратурной формулы
Для краткости обозначим
I [f ] = (b − a)s∑
k=0
wk f (xk), E [f ] =
∫baf (x)dx − I [f ]
Функционал E [f ] назовем остаточным членом квадратуры или просто,погрешностью квадратуры.Если квадратура имеет алгебраическую степень m, то E [Pm] ≡ 0.
Уткин П.С. Интегрирование 21 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член квадратурной формулы
Для определения E [f ] представим f (x) в виде формулы Тейлора синтегральным остаточным членом:
f (x) = f (a) + f ′(a)(x − a) + · · ·+ f (m)(a)
m!(x − a)m + Rm(x)
Rm(x) =
∫ xa
f (m+1)(t)
m!(x − t)mdt =
∫ba
f (m+1)(t)
m!(x − t)m+dt,
где2
z+ = max(z , 0) =
0, z 6 0z , z > 0
Так как f (x) − Rm(x) является многочленом степени m, онинтегрируется квадратурной формулой точно:
E [f ] = E [Rm]
2Условимся, что 00 = 0Уткин П.С. Интегрирование 22 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член квадратурной формулы
Остаточный член квадратурной формулы совпадает с погрешностьюинтегрирования остаточного члена формулы Тейлора.
Rm(x) =
∫ba
f (m+1)(t)
m!(x − t)m+dt,
E [f ] = E [Rm] =
∫baRm(x)dx − (b − a)
s∑k=0
wkRm(xk).
Подставим выражение для Rm(x) и поменяем последовательностьинтегрирования:
E [f ] =
∫ba
f (m+1)(t)
m!
[∫ba(x − t)m+dx − (b − a)
s∑k=0
wk(xk − t)m+
]︸ ︷︷ ︸
Q(t)
dt
Уткин П.С. Интегрирование 23 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член квадратурной формулы
Итак, остаточный член квадратуры выражен в виде
E [f ] =
∫ba
f (m+1)(t)
m!Q(t)dt,
где
Q(t) =
∫bt(x − t)mdx − (b − a)
∑k : xk>t
wk(xk − t)m
— функция, зависящая только от конкретного вида квадратурнойформулы и ее степени m. Заметим, что Q(t) достаточно вычислитьтолько для стандартного отрезка τ ∈ [−1, 1].
Q
(a + b
2+ τ
b − a
2
)=
(b − a
2
)m+1
q(τ)
q(τ) =
∫1τ
(ξ− τ)mdξ− 2∑
k : ξk>τ
wk(ξk − τ)m
Уткин П.С. Интегрирование 24 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член квадратурной формулы
Если Q(t) > 0 (< 0) на всем отрезке [a, b], то из теоремы обинтегральном среднем,
E [f ] =f (m+1)(ζ)
m!
∫baQ(t)dt =
f (m+1)(ζ)(b − a)m+2
2m+2m!
∫1−1
q(τ)dτ, ζ ∈ [a, b]
|E [f ]| 6Mm+1(b − a)m+2
2m+2m!
∣∣∣∣∫1−1
q(τ)dτ
∣∣∣∣Если Q(t) меняет знак, то неулучшаемая оценка выглядит так:
|E [f ]| 6Mm+1(b − a)m+2
2m+2m!
∫1−1
|q(τ)|dτ.
Уткин П.С. Интегрирование 25 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточные члены стандартных формул
Формула m q(τ) E [f ]
Левых прямоугольников 0 1− τ f ′(ζ)(b−a)2
2
Правых прямоугольников 0 −1− τ − f ′(ζ)(b−a)2
2
Средней точки 1 (1−|τ|)2
2f ′′(ζ)(b−a)3
24
Трапеций 1 τ2−12 − f ′′(ζ)(b−a)3
12
Симпсона 3 (|τ|−1)3(1+3|τ|)12 − f IV (ζ)(b−a)5
2880
Правило «3/8» 3
9τ4−1
36 , |τ| 6 1/3(|τ|−1)3|τ|
4 , |τ| > 1/3− f IV (ζ)(b−a)5
6480
Уткин П.С. Интегрирование 26 / 45
Численное интегрирование Погрешность квадратурных формул
Остаточный член составной формулы
Для элементарных квадратурных формул степени m остаточный членимеет вид
E [f ] =f (m+1)(ζ)(b − a)m+2
Cлибо |E [f ]| 6
Mm+1(b − a)m+2
C.
Найдем соответствующее выражение для составной квадратурнойформулы с равномерным шагом ∆xi = h = const:
E [f ] =N∑i=1
f (m+1)(ζi )hm+2
C= (b − a)
f (m+1)(ζ)hm+1
C, ζ ∈ [a, b].
|E [f ]| 6N∑i=1
Mm+1∆xm+2i
C6
Mm+1 max∆xm+1i
C
N∑i=1
∆xi = (b−a)Mm+1h
m+1
C.
Последнее выражение верно и для неравномерной сетки, еслиh ≡ max∆xi — максимальный шаг сетки.
Уткин П.С. Интегрирование 27 / 45
Численное интегрирование Погрешность квадратурных формул
Порядок сходимости
Теперь ясна связь между алгебраическим порядком квадратурнойформулы и порядком сходимости:
Утверждение 5.3.Для составной квадратурной формулы алгебраической степени m ееостаточный член стремится к нулю как O(hm+1), где h — шаг сеткисоставной формулы. Иначе говоря, формула имеет порядоксходимости m + 1.
Недостаточная гладкостьЕсли m + 1-я производная функции f (x) не ограничена, остаточныйчлен квадратурной формулы m-ой степени может стремиться к нулюмедленнее, чем O(hm+1). Остаточный член может быть найден тем жеспособом, но с заменой m на m ′ < m, где m ′ + 1 — числоограниченных производных у f (x).
Уткин П.С. Интегрирование 28 / 45
Численное интегрирование Погрешность квадратурных формул
Оценки остаточного члена при недостаточной гладкости
Построим для формулы Симпсона (m = 3) оценку остаточного членадля функции f (x) только с тремя (m ′ = 2) ограниченнымипроизводными. Построим функцию q(τ):
q(τ) =
∫1τ
(ξ− τ)m′dξ− 2
∑xik>τ
wk(ξk − τ)m ′
=
=(1− τ)3
3−
13(1− τ)2 −
43(−τ)2+ = −
τ(1− τ)2
3−
13(|τ|− τ)2 =
= −τ(|τ|− 1)2
3
Функция q(τ) не знакопостоянна на [−1, 1], и, следовательно,
|E [f ]| ≤ M3(b − a)4
24 · 2!
∫1−1
|q(τ)|dτ =M3(b − a)4
576
Уткин П.С. Интегрирование 29 / 45
Численное интегрирование Погрешность квадратурных формул
Оценки остаточного члена при недостаточной гладкости
Построим для формулы Симпсона (m = 3) оценку остаточного членадля функции f (x) только с двумя (m ′ = 1) ограниченнымипроизводными. Построим функцию q(τ):
q(τ) =
∫1τ
(ξ− τ)m′dξ− 2
∑xik>τ
wk(ξk − τ)m ′
=
=(1− τ)2
2−
13(1− τ) −
43(−τ)+ =
(1− 3τ)(1− τ)6
−23(|τ|− τ) =
=16+τ2
2− 2
|τ|
3
Функция q(τ) не знакопостоянна на [−1, 1], и, следовательно,
|E [f ]| ≤ M2(b − a)3
23 · 1!
∫1−1
|q(τ)|dτ =M2(b − a)3
81
Уткин П.С. Интегрирование 30 / 45
Численное интегрирование Апостериорное определение погрешности
Апостериорные оценки
• Построенные формулы и оценки для остаточного членаинтегрирования являются априорными3. На практике оценитьмаксимум нужной производной бывает сложно или вообщеневозможно (например, функция f (x) — «черный ящик»).
• Широкое распространение получили способы оценить ошибкупрямо в процессе вычислений. Оценки такого типа называютсяапостериорными4.
• Пусть I ∗ — точное значение интеграла, а Ih — приближенноезначение, вычисленное на сетке с шагом h по квадратурнойформуле порядка p (то есть степени p − 1):
I ∗ = Ih + E [f ] = Ih + (b − a)f (p)(ζ)hp
C
3лат. a priori — «от предшествующего», оценки известны до вычислений.4лат. a posteriori — «из последующего»
Уткин П.С. Интегрирование 31 / 45
Численное интегрирование Апостериорное определение погрешности
Правило Рунге
Сделаем предположение, что в выражении
I ∗ = Ih + E [f ] = Ih + (b − a)f (p)(ζ)hp
C
точка ζ слабо зависит от h:
I ∗ = Ih + chp + o(hp),
где c — константа и не зависит от h. Величина
εh = chp
может считаться главным членом ошибки (то есть ошибкой сточностью до бесконечно малой по h поправки).
Уткин П.С. Интегрирование 32 / 45
Численное интегрирование Апостериорное определение погрешности
Правило Рунге
Вычислим интеграл несколько раз на серии сгущающихся сеток сшагами h, h/2, h/4, . . . :
I ∗ = Ih + εh + o(hp) = Ih + chp + o(hp)
I ∗ = Ih/2 + εh/2 + o(hp) = Ih/2 + 2−pchp + o(hp)
I ∗ = Ih/4 + εh/4 + o(hp) = Ih/4 + 2−2pchp + o(hp)
...
Заметим, что разность приближенных значений интегралов позволяетоценить разность главных членов соответствующих погрешностей:
∆h/2 ≡ Ih/2 − Ih = εh − εh/2 + o(hp) = (2p − 1)εh/2 + o(h).
Уткин П.С. Интегрирование 33 / 45
Численное интегрирование Апостериорное определение погрешности
Правило Рунге
Правило Рунге позволяет просто оценить главный член погрешностивычислении интеграла на мелкой сетке:
εh/2 =Ih/2 − Ih
2p − 1+ o(h).
Использование правила Рунге требует осторожности. Требуетсяконтролировать что фактический порядок сходимости численногометода соответствует номинальному p. Фактический порядок p∗
может в силу некоторых обстоятельств (например, недостаточнаягладкость f (x)) быть меньше номинального порядка сходимости p.Простейший способ контроля — следить за выполнением соотношения
p∗ ≡ log2∆h
∆h/2≈ p.
Уткин П.С. Интегрирование 34 / 45
Численное интегрирование Апостериорное определение погрешности
Использование правила Рунге, пример 1
Интегрирование методом Симпсона p = 4 функции 1√xна [1, 9].
N Ih ∆h = Ih − I2h εh = ∆h/(24 − 1) p∗
40 4.0000010223489 ? ? ?
80 4.0000000647720 −9.57577 · 10−7 −6.38385 · 10−8 ?
160 4.0000000040624 −6.07096 · 10−8 −4.04731 · 10−9 3.98320 4.0000000002541 −3.80827 · 10−9 −2.53885 · 10−10 3.99640 4.0000000000159 −2.38246 · 10−10 −1.58831 · 10−11 4.001280 4.0000000000010 −1.48832 · 10−11 −9.92214 · 10−13 4.002560 4.0000000000001 −8.91731 · 10−13 −5.94487 · 10−14 4.065120 4.0000000000000 −1.16351 · 10−13 −7.75676 · 10−15 2.9410240 3.9999999999999 −1.28342 · 10−13 −8.55612 · 10−15 −0.14
Точное значение I =
∫91
dx√x= 4.
Уткин П.С. Интегрирование 35 / 45
Численное интегрирование Апостериорное определение погрешности
Использование правила Рунге, пример 2
Интегрирование методом Симпсона p = 4 функции 3−√x на [0, 9].
N Ih ∆h = I2h − Ih εh = ∆h/(24 − 1) p∗
40 9.0030633904588 ? ? ?
80 9.0010830724831 −1.98032 · 10−3 −1.32021 · 10−4 ?
160 9.0003829239736 −7.00149 · 10−4 −4.66766 · 10−5 1.50320 9.0001353840708 −2.47540 · 10−4 −1.65027 · 10−5 1.50640 9.0000478654974 −8.75186 · 10−5 −5.83457 · 10−6 1.501280 9.0000169230090 −3.09425 · 10−5 −2.06283 · 10−6 1.502560 9.0000059831870 −1.09398 · 10−5 −7.29321 · 10−7 1.505120 9.0000021153757 −3.86781 · 10−6 −2.57854 · 10−7 1.5010240 9.0000007478989 −1.36748 · 10−6 −9.11651 · 10−8 1.50
Точное значение I =
∫903−√xdx = 9.
Уткин П.С. Интегрирование 36 / 45
Численное интегрирование Формулы наивысшей степени
Степень интерполяционных формул
Утверждение 5.4Формулы интерполяционного типа, построенные по s + 1 узлугарантированно будут иметь алгебраическую степень не менее s.
Доказательство.На сетке из s + 1 узла интерполяция точно восстанавливаетмногочлены до степени s включительно. Следовательно, иквадратурная формула будет интегрировать их точно.
Формулы Ньютона-Котеса с нечетным количеством узлов имеютстепень на 1 выше, чем гарантируется утверждением 5.4:
Формула Количество узлов, s + 1 Алгебраическая степень, m
Трапеций 2 1 = s
Симсона 3 3 > s
Правило «3/8» 4 3 = sУткин П.С. Интегрирование 37 / 45
Численное интегрирование Формулы наивысшей степени
Условия повышения степени
Пусть I [f ] =s∑
k=0
wk f (xk) — интерполяционная квадратурная формула.
Рассмотрим произвольный многочлен P(x), степени выше s. Разделимего с остатком на многочлен ω(x) =
∏sk=0(x − xk):
P(x) = V (x)ω(x) + Z (x).
Многочлен Z (x) имеет степень не выше s, и, следовательно,интегрируется точно E [Z ] = 0. Тогда ошибка интегрирования P(x)составляет
E [P] = E [Vω] =
∫baV (x)ω(x)dx − I [Vω] =
∫baV (x)ω(x)dx
Квадратурная формула для V (x)ω(x) дает 0 поскольку эта функцияобнуляется во всех узлах xk .
Уткин П.С. Интегрирование 38 / 45
Численное интегрирование Формулы наивысшей степени
Условия повышения степени
Для многочленаP(x) = V (x)ω(x) + Z (x)
ошибка квадратурной формулы составляет
E [P] =
∫baV (x)ω(x)dx .
Например, для формулы Симпсона последний интеграл обнуляется,если V (x) = V — константа∫b
aV (x)(x − a)(x − b)
(x −
a + b
2
)dx =
= V
∫ba(x − a)(x − b)
(x −
a + b
2
)dx = 0
Стало быть, формула Симпсона будет точной для всех многочленов,которые при делении на ω(x) дают в частном константу, то есть длявсех многочленов степени s + 1.
Уткин П.С. Интегрирование 39 / 45
Численное интегрирование Формулы наивысшей степени
Условия повышения степени
Причина повышения порядка формулы Симпсона — в удачномрасположении узлов квадратуры.
ЗадачаРазместить узлы xk , k = 0, . . . , s квадратуры так, чтобы интеграл∫b
aV (x)ω(x)dx
обнулялся для произвольного многочлена V (x) как можно большейстепени d .
Если удастся так расположить узлы квадратуры, то она будет иметьалгебраический порядок m = s + d + 1.
Уткин П.С. Интегрирование 40 / 45
Численное интегрирование Формулы наивысшей степени
Ортогональные многочлены
Построим такой многочлен ω(x), который был бы ортогоналенмногочленам 1, x , . . . , x s в смысле
(x r ,ω(x)) ≡∫bax rω(x)dx = 0, r = 0, . . . , s
Для построения такого многочлена ω(x) можно использовать процессортогонализации Грама-Шмидта для набора функций 1, x , . . . , x s+1.В результате ортогонализции из многочлена x s+1 получится многочленω(x) = x s+1 + . . . , ортогональный всем функциям 1, x , . . . , x s , аследовательно, и любому многочлену V (x) степени degV 6 s.
Уткин П.С. Интегрирование 41 / 45
Численное интегрирование Формулы наивысшей степени
Многочлены Лежандра
-1.0 -0.5 0.5 1.0x
-1.0
-0.5
0.5
1.0
LHxL
L0HxL
L1HxL
L2HxL
L3HxL
L4HxL
Многочлены Лежандра Lm(x) ортогональны∫1−1
Lm(x)Ln(x)dx =2
n + 1δnm.
Каждый многочлен Lm(x) имеет m различных действительных корнейна отрезке [−1, 1].
Уткин П.С. Интегрирование 42 / 45
Численное интегрирование Формулы наивысшей степени
Многочлены Лежандра
Построенный многочлен ω(x) = Cs+1Ls+1(2x−a−b
b−a
), где Cs+1 —
некоторый несущественный нормировочный множитель.Использование в качестве узлов квадратурной формулы корнеймногочлена Лежандра позволяет добиться обнуления всех интегралов∫b
aV (x)ω(x)dx
при degV (x) 6 s, то есть повысить степень квадратуры до 2s + 1, тоесть почти вдвое по сравнению с формулами с произвольным выборомузлов.Для degV (x) = s + 1 никакой выбор узлов не обеспечит обнуленияинтеграла (достаточно взять V (x) = ω(x)). Стало быть,максимальная степень квадратуры с s + 1 узлом — 2s + 1.
Уткин П.С. Интегрирование 43 / 45
Численное интегрирование Формулы наивысшей степени
Формулы Гаусса
Формулами Гаусса называются интерполяционные квадратурныеформулы, имеющие максимальную алгебраическую степень дляданного числа узлов. Формула Гаусса с K узлами имеет степень2K − 1 и порядок 2K . Формулы Гаусса обычно приводят длястандартного отрезка [−1, 1]:
Число узлов, K Узлы xk Веса wk E [f ]
1 0 1 f ′′(ζ)(b−a)3
24
2 ± 1√3
1/2 f IV (ζ)(b−a)5
4320
30 4/9 f VI (ζ)(b−a)7
2016000±√
155 5/18
Простейшая формула Гаусса совпадает с формулой средней точки.Узлы формул Гаусса не содержат крайних точек отрезкаинтегрирования.
Уткин П.С. Интегрирование 44 / 45
Численное интегрирование Заключение
Список использованных источников
1 Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы.3-е издание, Москва. БИНОМ, Лаборатория знаний, 2004, 636 с.
2 Н.Н. Калиткин, Е.А. Альшина. Численные методы: в 2 кн. Кн. 1Численный анализ. Москва, Издательский центр «Академия»,2013, 304 с.
3 И.П. Мысовских. Интерполяционные кубатурные формулы.Москва. Наука, Главная редакция физико-математическойлитературы, 1981, 336 с.
Уткин П.С. Интегрирование 45 / 45