+ All Categories
Home > Documents > к.ф.-м.н.УткинПавелСергеевич e-mail:[email protected],pavel … ·...

к.ф.-м.н.УткинПавелСергеевич e-mail:[email protected],pavel … ·...

Date post: 28-Jul-2020
Category:
Upload: others
View: 9 times
Download: 0 times
Share this document with a friend
46
Численное интегрирование к.ф.-м.н. Уткин Павел Сергеевич 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]
Transcript
Page 1: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование

к.ф.-м.н. Уткин Павел Сергеевич 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]

Page 2: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Введение

Постановка задачи

Для заданной функции f (x) вычислить значение определенногоинтеграла

I =

∫baf (x)dx .

Желательно, чтобы метод численного интегрирования обладалследующими свойствами:• Универсальность. Функция f (x) может быть задана в виде«черного ящика», как способ вычисления f (x) по данному x .

• Экономичность. Количество вычислений функции f (x) повозможности должно быть сведено к минимуму.

• Хорошая обусловленность. Неустранимые погрешности ∆f взначениях f (x) не должны приводить к значительной итоговойошибке ∆I .

Уткин П.С. Интегрирование 2 / 45

Page 3: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Введение

Приложения

Численное интегрирование может применяться для• интегрирования функций, известных только в некоторых точках,например, полученных в результате измерений;

• интегрирования сложных выражений, не имеющих элементарныхпервообразных, либо имеющих слишком громоздкие выражениядля них;

• построения методов численного решения уравнений вобыкновенных и частных производных (методы конечныхэлементов, интегро-интерполяционные методы).

Уткин П.С. Интегрирование 3 / 45

Page 4: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Простейшие методы

Методы, основанные на определении интеграла

Вспомним, как определяется определенный интеграл Римана.Рассмотрим на отрезке интегрирования [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

Page 5: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Простейшие методы

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

Возьмем в качестве приближенного значения интеграла значениенекоторой интегральной суммы:∫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

Page 6: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Простейшие методы

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

∫baf (x)dx ≈

N∑i=1

f

(xi−1 + xi

2

)∆xi .

Уткин П.С. Интегрирование 6 / 45

Page 7: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Простейшие методы

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

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

∫baf (x)dx ≈

N∑i=1

f (xi−1)∆xi .

∫baf (x)dx ≈

N∑i=1

f (xi )∆xi .

Уткин П.С. Интегрирование 7 / 45

Page 8: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Простейшие методы

Составные квадратурные формулы

Построенные формулы прямоугольников являются составнымиквадратурными формулами. Квадратурная формула называетсясоставной, если является результатом применения некоторойэлементарной квадратурной формулы к каждому интервалу [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

Page 9: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Приближение подынтегральной функции

Пусть подынтегральная функция f (x) хорошо приближается некоторойпросто интегрируемой функцией P(x). Тогда∫b

af (x)dx ≈

∫baP(x)dx .

В качестве функции P(x) могут выступать• многочлены;• тригонометрические многочлены;• экспоненциальные многочлены и тому подобные.

Остановимся на случае, когда P(x) — многочлен. Тогда задачаприближения f (x) с помощью функции P(x) может быть решена,например, с помощью алгебраической интерполяции.

Уткин П.С. Интегрирование 9 / 45

Page 10: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Интерполяционные квадратурные формулы

Введем на отрезке [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

Page 11: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Интерполяционные квадратурные формулы

∫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

Уткин П.С. Интегрирование 11 / 45

Page 12: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Интерполяционные квадратурные формулы

Вид квадратурной формулы∫baf (x)dx ≈ (b − a)

s∑k=0

wk f (xk)

является универсальным, и мог быть написан из общих соображенийлинейности квадратурной формулы по значениям подынтегральнойфункции (по аналогии с линейностью самого интеграла).Величины wk называются весами квадратурной формулы. Веса независят от конкретного отрезка интегрирования, и могут бытьвычислены на некотором «стандартном» отрезке. Обычно используютотрезок [0, 1] или [−1, 1].Соответственно xk называются узлами квадратурной формулы. Онитакже обычно приводятся для стандартного отрезка, а дляконкретного отрезка [a, b] они получаются линейным преобразованием.

Уткин П.С. Интегрирование 12 / 45

Page 13: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Формулы Ньютона-Котеса

Изучим случай, когда Ω — равномерная сетка:

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

Page 14: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Формула трапеций

Для одного отрезка:∫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

Page 15: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Формула Симпсона

Для одного отрезка:∫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

Page 16: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Интерполяционные квадратурные формулы

Формула «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

Page 17: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Что точнее?

Попробуйте угадать, какая из двух квадратурных формул точнее:

O x

y

O x

y

Метод средней точки

Ошибка ∼ 1%

Метод трапеций

Ошибка ∼ 2%

Уткин П.С. Интегрирование 17 / 45

Page 18: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Что точнее?

Попробуйте угадать, какая из двух квадратурных формул точнее:

O x

y

O x

y

Метод средней точкиОшибка ∼ 1%

Метод трапецийОшибка ∼ 2%

Уткин П.С. Интегрирование 17 / 45

Page 19: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Степень квадратурной формулы

Будем говорить, что квадратурная формула имеет алгебраическуюстепень точности 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

Page 20: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Условия заданной степени

Утверждение 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

Page 21: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Условия заданной степени

Условия из утверждения 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

Page 22: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член квадратурной формулы

Для краткости обозначим

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

Page 23: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член квадратурной формулы

Для определения 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

Page 24: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член квадратурной формулы

Остаточный член квадратурной формулы совпадает с погрешностьюинтегрирования остаточного члена формулы Тейлора.

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

Page 25: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член квадратурной формулы

Итак, остаточный член квадратуры выражен в виде

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

Page 26: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член квадратурной формулы

Если 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

Page 27: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточные члены стандартных формул

Формула 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

Page 28: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Остаточный член составной формулы

Для элементарных квадратурных формул степени 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

Page 29: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Порядок сходимости

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

Утверждение 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

Page 30: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Оценки остаточного члена при недостаточной гладкости

Построим для формулы Симпсона (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

Page 31: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Погрешность квадратурных формул

Оценки остаточного члена при недостаточной гладкости

Построим для формулы Симпсона (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

Page 32: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Апостериорные оценки

• Построенные формулы и оценки для остаточного членаинтегрирования являются априорными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

Page 33: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Правило Рунге

Сделаем предположение, что в выражении

I ∗ = Ih + E [f ] = Ih + (b − a)f (p)(ζ)hp

C

точка ζ слабо зависит от h:

I ∗ = Ih + chp + o(hp),

где c — константа и не зависит от h. Величина

εh = chp

может считаться главным членом ошибки (то есть ошибкой сточностью до бесконечно малой по h поправки).

Уткин П.С. Интегрирование 32 / 45

Page 34: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Правило Рунге

Вычислим интеграл несколько раз на серии сгущающихся сеток сшагами 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

Page 35: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Правило Рунге

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

εh/2 =Ih/2 − Ih

2p − 1+ o(h).

Использование правила Рунге требует осторожности. Требуетсяконтролировать что фактический порядок сходимости численногометода соответствует номинальному p. Фактический порядок p∗

может в силу некоторых обстоятельств (например, недостаточнаягладкость f (x)) быть меньше номинального порядка сходимости p.Простейший способ контроля — следить за выполнением соотношения

p∗ ≡ log2∆h

∆h/2≈ p.

Уткин П.С. Интегрирование 34 / 45

Page 36: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Использование правила Рунге, пример 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

Page 37: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Апостериорное определение погрешности

Использование правила Рунге, пример 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

Page 38: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Степень интерполяционных формул

Утверждение 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

Page 39: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Условия повышения степени

Пусть 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

Page 40: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Условия повышения степени

Для многочлена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

Page 41: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Условия повышения степени

Причина повышения порядка формулы Симпсона — в удачномрасположении узлов квадратуры.

ЗадачаРазместить узлы xk , k = 0, . . . , s квадратуры так, чтобы интеграл∫b

aV (x)ω(x)dx

обнулялся для произвольного многочлена V (x) как можно большейстепени d .

Если удастся так расположить узлы квадратуры, то она будет иметьалгебраический порядок m = s + d + 1.

Уткин П.С. Интегрирование 40 / 45

Page 42: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Ортогональные многочлены

Построим такой многочлен ω(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

Page 43: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Многочлены Лежандра

-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

Page 44: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Многочлены Лежандра

Построенный многочлен ω(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

Page 45: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Формулы наивысшей степени

Формулы Гаусса

Формулами Гаусса называются интерполяционные квадратурныеформулы, имеющие максимальную алгебраическую степень дляданного числа узлов. Формула Гаусса с 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

Page 46: к.ф.-м.н.УткинПавелСергеевич e-mail:utkin@icad.org.ru,pavel … · 2014-10-04 · Численное интегрирование Простейшие методы

Численное интегрирование Заключение

Список использованных источников

1 Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы.3-е издание, Москва. БИНОМ, Лаборатория знаний, 2004, 636 с.

2 Н.Н. Калиткин, Е.А. Альшина. Численные методы: в 2 кн. Кн. 1Численный анализ. Москва, Издательский центр «Академия»,2013, 304 с.

3 И.П. Мысовских. Интерполяционные кубатурные формулы.Москва. Наука, Главная редакция физико-математическойлитературы, 1981, 336 с.

Уткин П.С. Интегрирование 45 / 45


Recommended