Новые многофункциональные материалы и нанотехнологии

114
Чаще всего приближенная функция выбирается в виде полинома, из- за простоты его интегрирования. Первообразная полинома степени m есть полином степени (m + 1).
Рассмотрим ситуацию, когда подынтегральная функция задана набором значений y
i
для произвольных узлов x
i
, принадлежащих интервалу [a, b], причем i = 0, 1, …, n. Оба вышеупомянутых варианта сводятся к данной постановке задачи.
Согласно изложенному в главе 2, для произвольного набора (n+1) точек x
i
, y
i
(i = 0, 1, …, n) существует единственный интерполяционный полином Лагранжа (2.3). Этот полином L
n
(x) подставим в (8.1) на место подынтегральной функции. При этом мы получаем приближенное значение определенного интеграла:
I
(f, a, b)
≈
∫
b
a
n
dx
x
L
)
(

(8.2)
Применим теорему Ньютона–Лейбница и представим результат в виде
I(f, a, b)
≈
∑
=
n
i
i
i
y
A
0
,

(8.3) где A
i

(i = 0, 1, …, n) – постоянные коэффициенты, полученные в результате интегрирования полинома Лагранжа (2.3) и подстановки пределов:
A
i

=
∫
b
a
i
dx
x
)
(
l
,

(8.4) где функции l
i

(i = 0, 1, …, n) представляются произведениями (2.7).
Заметим, что функции l
i

(i = 0, 1, …, n) представляют собой полиномы
n
-й степени, полностью определяемые положением узлов x
i
(i = 0, 1, …,
n
) на интервале интегрирования.
Следовательно, для приближенного вычисления определенного интеграла требуется вычислить сумму (8.3) произведений значений функции y
i
в узлах x
i
на коэффициенты A
i
, определяемые формулой

115
(8.4). Сумма, стоящая в правой части равенства (8.3), называется квадратурной, а выражение (8.3) называется приближенной
квадратурной формулой.
Таким образом, основная трудность численного интегрирования сводится к нахождению коэффициентов A
i

для квадратурной формулы
(8.3). Заметим, что значения коэффициентов A
i

(i = 0, 1, …, n), вычисляемых согласно (8.4) с помощью выражений (2.7), не зависят от вида интегрируемой функции f(x). Видно, что произведения (2.7) зависят только от расположения узлов x
i
(i = 0, 1, …, n) внутри интервала интегрирования
[a,
b
].
Из-за единственности интерполяционного полинома Лагранжа L
n
(x), проведенного через заданные точки интегрируемой функции (x
i
,
y
i
), значение определенного интеграла, рассчитанного по квадратурной формуле
(8.3), определяется только расположением узлов x
i
(i = 0, 1, …, n).
Проводить практическое вычисление коэффициентов
A
i
квадратурной формулы (8.3) путем интегрирования полиномов n-й степени l
i
, которые представляются произведениями (2.7), довольно затруднительно. Гораздо более простой способ основан на том, что для подынтегральных функций f(x) = x
k
, где k = 0, 1, …, n, квадратурная формула (8.3) является точной. Иначе говоря, равенства
∫
b
a
k
dx
x
=
∑
=
n
i
i
i
y
A
0
, k = 0, 1, …, n
(8.5) справедливы для любого расположения узлов x
i
(i = 0, 1, …, n) на интервале интегрирования [a, b].
С другой стороны, определенные интегралы в правых частях уравнений (8.5) вычисляются по формуле Ньютона–Лейбница:
∫
b
a
k
dx
x
=
1 1
1
+
−
+
+
k
a
b
k
k
, k = 0, 1, …, n .
(8.6)
Приравнивая правые части соответствующих уравнений (8.5) и (8.6), получим n + 1 уравнений для вычисления n + 1 искомых коэффициентов

116
A
i
(i = 0, 1, …, n). Систему полученных уравнений можно записать в виде:
I
0

=
∑
=
n
i
i
A
0
, I
1

=
∑
=
n
i
i
i
x
A
0
, … , I
n

=
∑
=
n
i
n
i
i
x
A
0
,
(8.7) где числовые значения правых частей уравнений (8.6) обозначены I
k
(k =
0, 1, …, n).
Система уравнений (8.7) является линейной относительно неизвестных величин коэффициентов A
i
(i = 0, 1, …, n). Способы решения таких систем описаны в главе 4. Таким образом, вычисляются искомые коэффициенты A
i
(i = 0, 1, …, n) и квадратурная формула (8.3) готова к применению. Следует отметить, что, хотя при получении квадратурной формулы использовался интерполяционный полином
Лагранжа (2.3), описанный способ численного интегрирования не требует построения этого полинома в явном виде.
Пример 1. Пусть интервал интегрирования определяется границами
a
= 0 и b = 1. Известны значения подынтегральной функции в узлах
x
0
=1/4, x
1
= 1/2, x
2
= 3/4. Подсчет величин I
k
(k = 0, 1, 2) по формулам
(8.6) дает числовые значения I
0
= 1, I
1
= 1/2, I
2
= 1/3. Система (8.7) примет вид:
A
0
+ A
1
+ A
2
= 1,
A
0
/ 4 + A
1
/ 2 + 3 A
2
/ 4 = 1/2,
A
0
/ 16 + A
1
/ 4 + 9 A
2
/ 16 = 1/3.
Решением последней системы являются следующие значения коэффициентов A
0
= 2/3, A
1
=
−1/3, A
2
= 2/3. Теперь вычисление приближенного значения определенного интеграла может проводиться по следующей квадратурной формуле:
∫
1 0
)
( dx
x
f
≈ 2 y
0
/ 3
− y
1
/ 3 + 2 y
2
/ 3, где y
0
= f (x
0
= 1/4), y
1
= f (x
1
= 1/2), y
2
= f (x
2
= 3/4).

117
Примечание. Описанный в данном параграфе способ вычисления коэффициентов квадратурной формулы базируется на замене подынтегральной функции f (x) интерполирующим полиномом. Если через заданные точки x
i
, y
i
(i = 0, 1, …, n) провести другую интерполирующую функцию (не полином), то приближенное вычисление определенного интеграла (8.1) также можно проводить по квадратурной формуле (8.3). Но при этом коэффициенты A
i
нельзя вычислять ни по формуле (8.4), ни решением системы (8.7). Значения коэффициентов A
i
и способ их вычисления зависят от вида выбранной интерполирующей функции.

§ 8.2. Квадратурные формулы Ньютона–Котеса
В предыдущем параграфе было показано, что численные значения коэффициентов A
i
(i = 0, 1, …, n) определяются положением узлов x
i
(i=0, 1, …, n) на диапазоне интегрирования [a, b]. Выбирая узлы x
i
различными способами, мы будем получать разные наборы коэффициентов A
i
Предположим, что мы можем вычислить n + 1 значений подынтегральной функции f(x) в произвольных n + 1 точках диапазона интегрирования [a, b].
Рассмотрим в первую очередь равномерное распределение узлов x
i
по диапазону интегрирования [a, b]. Для этого разобьем диапазон [a, b] на n равных интервалов, каждый длиной
h =
(b – a) / n.

(8.8)
Узлами будут границы интервалов:
x
i
= a + i h
, (i = 0, 1, …, n).
(8.9)
Ясно, что x
0
= a , x
n
= b
Вычислим (или измерим) значения подынтегральной функции в выбранных узлах y
i
= f
(x
i
) , i = 0, 1, …, n. Теперь, как и в предыдущем параграфе, заменим подынтегральную функцию f(x) интерполяционным полиномом Лагранжа (2.3). Коэффициенты квадратурной формулы A
i
(i=0, 1, …, n) будем вычислять по формуле (8.4) с помощью (2.7). Тогда

118 выражения для коэффициентов A
i
могут быть представлены в следующем виде:
A
i
=
(b – a) H
i
, i = 0, 1, …, n,
(8.10) где введены новые обозначения:
H
i

=
)!
(
!
)
1
(
i
n
i
n
i
n
−
−
−
dq
i
q
n
q
q
q
n
∫
−
−
−
0
)
)...(
1
(
,
(8.11)
q = (x – a) / h
,
dq = dx / h

Параметры (8.11) называются коэффициентами Котеса. Нетрудно убедиться, что при любом конечном n выполняются равенства:
H
i

= H
n – i

и
∑
=
n
i
i
H
0
= 1.
Таким образом, для вычисления определенного интеграла (8.1) следует применять квадратурную формулу следующего вида:
∫
b
a
dx
x
f
)
(
≈ (b – a)
∑
=
n
i
i
i
y
H
0

(8.12)
Формулы (8.12), где коэффициенты H
i

определены равенствами
(8.11), называются квадратурными формулами Ньютона–Котеса n-го порядка.
Для практического применения целесообразно рассмотреть некоторые частные случаи.
Пусть n = 1. Это значит, что диапазон интегрирования [a, b] содержит только два узла: x
0
= a
и x
1
= b
Вычислим коэффициенты Котеса для этого случая. Подставляя в
(8.11) n = 1 и i = 0, 1 получим:

119
H
0
= – 1
dq
q
q
q
∫
−
1 0
)
1
(
= 1 / 2, H
1

= 1
dq
q
∫
1 0
= 1 / 2.
(8.13)
Подстановка коэффициентов (8.13) в (8.12) дает квадратурную формулу Ньютона–Котеса первого порядка:
I
(f, a, b) =
∫
1 0
)
(
x
x
dx
x
f
≈ h (y
0
+ y
1
) / 2 = (b – a) (y
0
+ y
1
) / 2 . (8.14)
В данном случае величина шага h равна длине всего диапазона интегрирования [a, b].
При таким приближении подынтегральная функция f (x) заменяется линейной (полиномом Лагранжа 1-й степени). Графически значение определенного интеграла I(f, a, b) изображается площадью криволинейной трапеции, ограниченной сверху (или снизу) графиком функции f (x) (см. пример на рис. 8.1).







Рис. 8.1. К выводу формулы трапеций.
Сплошная линия – график интегрируемой функции f (x), штрихпунктирная – график линейного приближения интегрируемой функции на интервале [x
0
, x
1
]. Крестиками заполнена область R, площадь которой равна погрешности квадратурной формулы (8.14).
x
1

x
0

x

y
1

y
0

f
(x)

y

h

R

120
Приближенное значение интеграла, выраженное квадратурной суммой (8.14), равно площади обычной трапеции с той же высотой h.
Площадь криволинейного сектора дает погрешность R вычисления определенного интеграла с помощью квадратурной формулы (8.14).
Когда данный способ вычисления определенного интеграла используют на практике, то обычно диапазон интегрирования [a, b] разбивают на N одинаковых поддиапазонов шириной
h =
(b – a) / N

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

(k = 0, 1,
…
, N) в эту сумму будет входить дважды, кроме двух крайних y
0
и y
N
. В результате интеграл I(f, a, b) выразится следующей суммой:
I
(f, a, b)
≈ h {
∑
−
=
1 0
N
k
k
y +
2 0
N
y
y
+
}.

(8.16)
Последняя формула называется формулой трапеций.
Теперь пусть n = 2. Это значит, что диапазон интегрирования [a, b] содержит три узла: x
0
= a
, x
1
= (a + b) / 2, x
2
= b
. В данном случае шаг интегрирования (8.8) равен половине диапазона интегрирования
h
= (b – a) / 2.
Коэффициенты Котеса в данном случае получаются вычислением интегралов (8.11) с подстановкой n = 2 для i = 0, 1, 2:
H
0
=
2 1
2 1
dq
q
q
∫
−
−
2 0
)
2
)(
1
(
=
6 1
,
H
1
= –
2 1
2 1
dq
q
q
∫
−
2 0
)
2
(
=
3 2
,
(8.17)
H
2
=
2 1
2 1
dq
q
q
∫
−
2 0
)
1
(
=
6 1

Искомый интеграл представится в виде:

121
I
(f, a, b) =
∫
2 0
)
(
x
x
dx
x
f
≈ h (y
0
+ y
1
) / 2 = (b – a) (
6 1
y
0
+
3 2
y
1
+
6 1
y
2
), (8.18) где (b – a) = 2h.
Формула (8.18) была получена заменой подынтегральной функции
f
(x) на квадратичную, проходящую через заданные точки (x
i
, y
i
) , i = 0,
1,
2.
Полученное значение искомого интеграла графически изображается площадью под параболой (см. рис. 8.2).





Рис. 8.2. К выводу формулы Симпсона.
Жирная сплошная линия – график интегрируемой функции f (x), штриховая – график квадратичного приближения L
2
(x) интегрируемой функции на интервале [x
0
, x
2
]. Заштрихованная область R выражает погрешность квадратурной формулы (8.18). Видно, что на интервале [x
0
,
x
1
] квадратичное приближение дает заниженный результат, на интервале [x
1
, x
2
] – завышенный.
На практике весь диапазон интегрирования [a, b] обычно предварительно разбивается на четное число N = 2M равных интервалов, которые попарно объединяются в M поддиапазонов.
Каждый поддиапазон содержит три узла: два крайних и один
R
y
0

y
1

y
2

y

x
0

x
1

x
2

h

h

y=L
2
(x)

y=f
(x)

+

-

M
0
M
1
M
2
x

122 внутренний. К каждому из M поддиапазонов применяется формула
(8.18), т.е. каждой тройке узлов соответствует своя интерполирующая парабола. Параболы стыкуются в совпадающих крайних узлах отдельных поддиапазонов. Номера всех таких узлов четные от i = 2 до i
= 2 M
− 2, поэтому при суммировании выражений (8.18) по поддиапазонам следует учесть, что слагаемые с этими номерами встретятся дважды.
В результате суммирования с учетом значений коэффициентов
Котеса (8.17) получается рабочая формула
I
(f, a, b)
≈
3
h
(y
0
+ 4 S
O
+ 2 S
E
+ y
2М
),
(8.19) где для краткости введены обозначения следующих сумм:
S
O
= y
1
+ y
3
+ … + y
2М
−
1
, S
E
= y
2
+ y
4
+ … + y
2М
−
2
(8.20)
Заметим, что сумма S
O
содержит M слагаемых, а сумма S
E
состоит из M
−1 слагаемых. Шаг интегрирования равен h=(b – a)/(2M)=(b – a)/N.
Формула (8.19) называется квадратурной формулой Симпсона.
При n = 3 аналогичным путем получается квадратурная формула с тремя узлами на интервале интегрирования:
I
(f, a, b) =
∫
3 0
)
(
x
x
dx
x
f
≈
8 3
h
(y
0
+ 3y
1
+ 3y
2
+ y
3
),
(8.21) где h = (b – a) / 3.
При решении прикладных задач диапазон интегрирования [a, b] предварительно разбивается на равные поддиапазоны, число которых кратно трем (N = 3M), и к каждому из поддиапазонов применяется формула (8.20). Тогда получится рабочая формула вида:
I
(f, a, b)
≈
8 3
h
(y
0
+ 3 S
1
+ 3 S
2
+ 2 S
3
+ y
N
).
(8.22)
Шаг интегрирования равен h = (b – a) / (3M) = (b – a) / N, а отдельные суммы выражены следующими формулами:

123
S
1
=
∑
−
=
+
1 0
1 3
M
k
k
y
, S
2
=
∑
−
=
+
1 0
2 3
M
k
k
y
, S
3
=
∑
−
=
1 1
3
M
k
k
y

(8.23)
Сумма (8.22) называется квадратурной формулой Ньютона или формулой «трех восьмых»
. В практике квадратурная формула
Ньютона употребляется гораздо реже, чем формулы Симпсона и трапеций.
Еще реже применяется квадратурная формула Боде, которая получается из общей формулы (8.12) подстановкой n = 4. Интервал интегрирования содержит 5 равноотстоящих узлов x
i
= a + i h
, i = 0, 1,
…, 4. Приближенное значение определенного интеграла выразится следующей суммой:
I
(f, a, b)
≈
45 2
h
(7y
0
+ 32y
1
+ 12y
2
+ 32y
3
+ 7y
4
),
(8.24) где h = (b – a) / 4.
Читателю предоставляется возможность самостоятельно записать рабочую формулу численного интегрирования при разбиении диапазона интегрирования [a, b] на N равных поддиапазонов.

§ 8.3. Квадратурная формула Гаусса
Все формулы численного интегрирования дают результат с некоторой погрешностью. Оценки погрешностей приводятся в следующем параграфе. Здесь заметим, что уменьшение погрешности расчета является одной из важнейших задач численных методов.
В предыдущем параграфе рассматривались различные квадратурные формулы для равномерного расположения узлов в интервале интегрирования [a, b]. Оказывается, при фиксированном количестве узлов можно добиться значительного уменьшения погрешности вычисления определенного интеграла, если отказаться от их равномерного расположения.
Для нахождения оптимального расположения узлов (которое обеспечивает уменьшение погрешности численного интегрирования) прежде всего сведем интервал интегрирования общего вида [a, b] к

124 стандартному интервалу [–1, 1] с помощью линейной замены переменной интегрирования:
x
= (a + b) / 2 + t (b – a) / 2.
(8.25)
Исходный интеграл (8.1) приобретет вид:
∫
b
a
dx
x
f
)
(
=
∫
−






−
+
+
1 1
2 2
dt
t
a
b
b
a
f
,
(8.26) где t – новая переменная интегрирования.
В соответствии с общей квадратурной формулой (8.3) запишем:
( )
∫
−
1 1
dt
t
f
≈
∑
=
n
i
i
i
t
f
A
1
)
(

(8.27)
Нашей целью является выбор таких n узлов t
i
внутри интервала интегрирования [–1, 1] и таких n коэффициентов A
i
(i = 1, …, n), чтобы формула (8.27) была точной для всех полиномов максимально большой степени. Мы имеем возможность варьировать 2n величин t
i
и A
i
(i = 1,
…, n), следовательно, эта максимальная степень полинома равна
N
= 2n – 1.
Далее требуется доказать следующую лемму.
Лемма. Для достижения поставленной цели необходимо и достаточно, чтобы формула (8.27) выполнялась точно для следующих
2n подынтегральных функций f(t) = 1, t, t
2
, …, t
2n – 1
Доказательство.
По условиям леммы выполняются следующие равенства для степенных функций f(t) = t
k
(k = 0, 1, …, 2n – 1):
∫
−
1 1
dt
t
k
=
∑
=
n
i
k
i
i
t
A
1
, k = 0, 1, …, 2n – 1.
(8.28)
Как было принято выше, подынтегральную функцию f(t) представим в виде полинома степени 2n – 1:

125
f
(t) =
∑
−
=
1 2
0
n
k
k
k
t
C
,

(8.29) где С
k
– коэффициенты полинома.
Проведем интегрирование полинома (8.29) по интервалу [–1, 1], используя условие (8.28) и изменяя порядок суммирования:
∫
−
1 1
)
(
dt
t
f
=
∑
−
=
1 2
0
n
k
k
C
∫
−
1 1
dt
t
k
=
∑
−
=
1 2
0
n
k
k
C
∑
=
n
i
k
i
i
t
A
1
=
=
∑
=
n
i
i
A
1
∑
−
=
1 2
0
n
k
k
i
k
t
C
=
∑
=
n
i
i
i
t
f
A
1
)
(

(8.30)
Последнее равенство обусловлено определением (8.29).
Мы получили, что равенство (8.27) является точным для полинома вида (8.29), при этом в ходе преобразований мы воспользовались условиями (8.28).
Лемма доказана.
Эта лемма обеспечивает справедливость следующего важного утверждения. Пусть нам удастся найти такие числа t
i
и A
i
(i = 1, …, n), обращающие формулу (8.27) в точную, используя в качестве подынтегральных степенные функции вида t
k
(k = 0, 1, …, 2n – 1). Тогда эти найденные числа t
i
и A
i
обеспечат точность формулы (8.27) и для любой подынтегральной функции, которая представляется в виде полинома степени 2n – 1.
Теперь в уравнениях (8.28) проинтегрируем аналитически левые части:
∫
−
1 1
dt
t
k
=






−
−
+
)
(
,
0
)
(
,
1 2
нечетно
k
четно
k
k

(8.31)
Сравнение (8.31) с правыми частями уравнений (8.28) дает нам систему 2n уравнений для 2n искомых неизвестных t
i
, A
i
(i = 1, …, n):

126
∑
=
n
i
i
A
1
= 2,
∑
=
n
i
i
i
t
A
1
= 0,
………………………………………………
(8.32)
∑
=
−
n
i
n
i
i
t
A
1 2
2
=
1 2
2
−
n
,
∑
=
−
n
i
n
i
i
t
A
1 1
2
= 0.
Система является нелинейной относительно неизвестных t
i
, A
i
(i = 1,
…, n). Для ее решения целесообразно воспользоваться свойствами полиномов Лежандра. Краткие сведения о полиномах Лежандра приведены в приложении 2.
Выберем подынтегральные функции f(x) в виде:
f
(x) = t
k
P
n
(t), k = 0, 1, …, n – 1,
(8.33) где P
n

(t) - полином Лежандра степени n.
Это значит, что все n функций (8.33) являются полиномами степени
≤ 2n – 1. Для таких подынтегральных функций, согласно доказанному выше, формула (8.27) является точной и коэффициенты A
i
(i = 1, …, n) удовлетворяют системе (8.32).
Применим квадратурную формулу (8.27) для функций (8.33):
∫
−
1 1
)
(
dt
t
P
t
n
k
=
∑
=
n
i
i
n
k
i
i
t
P
t
A
1
)
(
, k = 0, 1, …, n – 1. (8.34)
В силу свойства ортогональности полиномов Лежандра левые части всех уравнений (8.34) равны нулю. Следовательно, и правые части
(8.34) будут равны нулю
∑
=
n
i
i
n
k
i
i
t
P
t
A
1
)
(
= 0

(8.35) при произвольных значениях A
i
, если положить P
n

(t
i
) = 0.
Таким образом, оказывается, что в качестве узлов t
i
следует взять точки нулей полиномов Лежандра степени n. Точки нулей полиномов

127
Лежандра различных степеней рассчитаны с большой точностью и сведены в специальные таблицы. Характерным свойством узлов t
i
и коэффициентов A
i
в данном методе является симметрия их значений относительно центра таблицы t = 0.
Если числа t
i

известны, то система (8.32) становится линейной относительно величин A
i
. Для ее решения можно применить любой метод, описанный в главе 4. Численные значения узлов t
i
и коэффициентов A
i
с точностью до восьми знаков для различных n
приведены в таблице приложения 3.
Вычислив значения узлов t
i
и коэффициентов A
i
для выбранного n, можно вернуться к исходной переменной интегрирования x и записать квадратурную

Скачать 0.96 Mb.

Поделитесь с Вашими друзьями: