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

128
§ 8.4. Погрешности квадратурных формул
Погрешность любой формулы численного интегрирования определяется остаточным членом, равным разности точного значения определенного интеграла и квадратурной суммы вида (8.3):
R
=
∫
b
a
dx
x
f
)
(
−
∑
=
n
i
i
i
y
A
0

(8.39)
Так как точное значение интеграла неизвестно, то остаточный член приходится вычислять приближенно с помощью известных значений подынтегральной функции y
i
= f (x
i
).
За погрешность численного интегрирования с помощью квадратурной формулы можно принять приближенное абсолютное значение остаточного члена (8.39) или его оценку сверху.
Обратимся вначале к методу трапеций. Интервал интегрирования [a,
b
] разбивается на N равных поддиапазонов шириной h, и к каждому из них применяется формула (8.14). Остаточный член (8.39) для любого поддиапазона может быть представлен в виде:
R
=
−
12 3
h
f
′′(ξ),

(8.40) где
ξ − некоторая точка внутри данного поддиапазона. Рисунок 8.1 показывает, что в случае использования метода трапеций абсолютная величина остаточного члена (8.39) равна площади криволинейного сектора. По рис. 8.1 хорошо видно, что при f
′′(x) < 0 на интервале интегрирования (x
0
, x
1
) квадратурная формула (8.14) дает значение с недостатком и остаточный член (8.40) имеет положительное значение.
При f
′′(x) > 0 формула (8.14) дает значение определенного интеграла с избытком и величина остаточного члена R согласно (8.40) отрицательна.
Интеграл по всему диапазону [a, b] дается квадратурной формулой трапеций (8.16). Погрешность этой формулы определяется суммой остаточных членов (8.40) по всем N поддиапазонам. Оценка сверху абсолютной величины остаточного члена формулы (8.16) может быть записана в следующем виде:

129
R
= h
2 12
a
b
−
M
2
,

(8.41) где M
2
− максимум абсолютного значения производной 2-го порядка функции f (x) на интервале интегрирования [a, b]:
M
2
=
b
x
a
≤
≤
max
 f ′′(x) .
(8.42)
Если отсутствует аналитическое выражение второй производной
f
′′(x), то для оценки величины M
2
используется представление конечной разностью (1.3) или формула численного дифференцирования, приводимая в следующей главе.
Перейдем к оценке погрешности численного интегрирования по методу Симпсона. В этом методе диапазон интегрирования [a, b] разбит на четное число N = 2M равных интервалов шириной h, к каждой паре из которых применяется квадратурной формула (8.18). Остаточный член этой формулы представляется так:
R
=
−
90 5
h
f
(IV)
(
ξ),

(8.43) где
ξ − некоторая точка внутри данного поддиапазона.
Численное интегрирование методом Симпсона по всему диапазону
[a, b] дается квадратурной формулой (8.19). Погрешность этой формулы определяется суммой остаточных членов (8.43) по всем M= N/2 поддиапазонам. Приближенное значение погрешности можно записать в виде:
R
= h
4 180
a
b
−
M
4
,

(8.44) где M
4
− максимум абсолютного значения производной 4-го порядка
f
(IV)
(x) подынтегральной функции f (x) на интервале [a, b]:
M
4
=
b
x
a
≤
≤
max
)
(
)
(
x
f
IV

(8.45)

130
Численная оценка величины f
(IV)
(x) при отсутствии аналитического вида функции f (x) весьма затруднительна. Поэтому на практике оценку погрешности вычисления определенного интеграла методом Симпсона реализуют эмпирическим путем. Интеграл (8.1) вычисляют по формуле
(8.19) дважды: с малым шагом h
1
и с вдвое меньшим шагом h
2
= h
1
/ 2.
Соответствующие значения квадратурных сумм (8.19) обозначим I
1
и I
2
Тогда в качестве оценки искомого интеграла берется величина I
2
Если интегрируемая функция достаточно гладкая, то ее четвертая производная f
(IV)
(x) на диапазоне [a, b] изменяется медленно и ее можно приближенно заменить константой. Тогда оценка погрешности метода
Симпсона представима в виде:
R
≈ |I
2
– I
1
| / 15.

(8.46)
Расчет погрешностей для квадратурных формул Ньютона
−−−−
Котеса более высоких порядков приводится в специальных курсах численных методов, так как эти формулы применяются гораздо реже, чем формулы трапеций и Симпсона.
Погрешность значения определенного интеграла по квадратурной формуле Гаусса также определяется величиной остаточного члена этой формулы. Остаточный член квадратурной формулы Гаусса резко зависит от количества узлов n в отдельном интервале интегрирования
[a, b]. В общем виде остаточный член, определяющий погрешность, представляется следующим выражением:
R
n
=
)
1 2
(
)
!
2
(
)
!
(
)
(
3 4
1 2
+
−
+
n
n
n
a
b
n
f
(2n)
(
ξ),
(8.47) где
ξ ∈ [a, b]. Например, для n = 3 и для n = 6 получаем:
R
3
=
15750 1
7 2






− a
b
f
(6)
(
ξ) и R
6
=
50 6489844861 1
13 2






− a
b
f
(12)
(
ξ).
Однако необходимость численных оценок производных высоких порядков делает формулу (8.47) мало пригодной для практики. Оценки погрешности приходится выполнять эмпирически, постепенно удваивая количество поддиапазонов интегрирования.

131
Пример 2. Для сравнения точности различных квадратурных формул рассмотрим в качестве примера вычисление значения функции erf(x=1).
Эта функция задана интегралом, где аргумент является верхним пределом:
erf
(x) =
π
2
∫
−
x
t
0 2
)
exp(
dt

(8.48)
Вычисление erf(x = 1) должно проводиться численно, так как для первообразной функции exp(
−t
2
) не существует аналитического выражения.
Расчет значения определенного интеграла
π
2
∫
−
1 0
2
)
exp( t
dt
проведем тремя способами: по 2-точечной формуле трапеций (8.14), по
3-точечной формуле Симпсона (8.18) и по 2-точечной формуле Гаусса
(8.36). При этом диапазон интегрирования не будем разбивать на поддиапазоны. Получим:
I
T
=
π
1
(1 + e
−1
) = 0,7717433… ,
I
S
=
π
3 1
(1 + 4 e
−1/4
+ e
−1
) = 0,8431028… ,
I
G
=
π
1






















+
−
+














−
−
2 2
3 1
1 4
1
exp
3 1
1 4
1
exp
= 0,8424419… . где I
T
, I
S
и I
G
− значения определенного интеграла, вычисленные методами трапеций, Симпсона и Гаусса соответственно.
Точный расчет значения функции erf(x = 1) с точностью до 10
–8
проводится с помощью специального ряда и дает 0,84270079… .
Сравнение результатов показывает, что значение, вычисленное по квадратурной формуле трапеций, значительно отличается от точного результата. Квадратурная формула Гаусса дает меньшую погрешность по сравнению с формулой Симпсона.
Конечно, уменьшить погрешность любого метода можно, разбивая диапазон [a, b] интегрирования на поддиапазоны и используя формулы

132
(8.16), (8.19) и (8.38). Однако различие в точности методов остается тем же.
Вернемся к задаче вычисления функции erf(x = 1). В таблице 8.1 представлены результаты численного интегрирования методом трапеций (8.16) по диапазону [0, 1] для разного числа интервалов.

Таблица 8.1
Число интервалов интегрирования
Значение квадратурной суммы
2 0,82526295 4
0,83836777 8
0,84161922 16 0,84243051 32 0,84263323 64 0,84268390 128 0,84269657
Видно, что при увеличении количества поддиапазонов метод трапеций демонстрирует весьма медленную сходимость к точному значению определенного интеграла. Поэтому метод трапеций обычно используют для быстрых приближенных оценок, которые содержат малое количество правильных значащих цифр.
Точность квадратурной формулы Гаусса выше по сравнению с формулой
Симпсона при заданном разбиении диапазона интегрирования. Это позволяет, применяя метод Гаусса, при заданной погрешности вычисления определенных интегралов разбивать диапазон интегрирования на меньшее число подинтервалов, что приводит к значительному сокращению времени расчета. Опыт показывает, что количество поддиапазонов при использовании квадратурной формулы
Гаусса может быть взято, грубо говоря, на порядок меньше, чем при применении метода Симпсона.
Примечание. В данном пособии не рассматриваются методы вычисления многомерных и несобственных интегралов. По этим вопросам авторы отсылают читателей к более подробным курсам численных методов, например к [3, 9].

133
Глава 9. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ
§ 9.1. Дифференцирование интерполяционных полиномов
Вычисление производных заданных функций проводится обычно по хорошо известным правилам и формулам дифференцирования. При этом получается аналитическое выражение для производной функции.
Подстановка в полученное выражение любого аргумента, принадлежащего области допустимых значений, позволяет вычислить искомое значение производной функции.
Численные (неаналитические) методы вычисления значений производной функции привлекаются в следующих случаях, встречающихся на практике.
Во-первых, когда исходная функция задана лишь таблицей значений в узлах. Промежуточные значения данной функции (для аргументов, не совпадающих с узлами) вычисляются методами интерполяции или аппроксимации.
Во-вторых, когда исходная функция определяется некоторым алгоритмом, который позволяет получить значение y = f(x) для любого аргумента x, принадлежащего области определения функции, но аналитическое выражение (формула) для функции f(x) отсутствует.
Следовательно, стандартные формулы дифференцирования неприменимы.
В-третьих, зависимость y = f(x) может выражаться сложной формулой с использованием специальных функций. Нахождение аналитического представления производной сопряжено с изрядными трудностями, а формула f
′(x) приобретает весьма громоздкий вид, мало удобный для вычислений. В частности, выражение может представляться бесконечным рядом неэлементарных функций, значения которых в принципе вычисляются приближенно. Во многих подобных случаях численные методы дифференцирования имеют преимущество – они позволяют получить значение производной с достаточной точностью за меньшее время по сравнению с расчетом по полученной формуле f
′(x).
Прежде чем переходить к конкретным методам, следует обратить внимание на принципиальную трудность численного дифференцирования. В курсе математического анализа доказано, что для сколь угодно близких непрерывных функций расстояние между их производными может быть неограниченно велико. Это означает, что

134 задача дифференцирования на пространстве непрерывных функций в принципе некорректна.
Вышеуказанное свойство функций обусловливает ограниченную точность всех формул численного дифференцирования.
Целью методов численного дифференцирования является расчет значения производной f
′(x) функции y = f (x) в определенной точке x.
Обычные методы численного дифференцирования базируются на имеющемся алгоритме вычисления значений исходной функции y = f (x) для заданного аргумента x. В первом из перечисленных случаев функция f(x) заменяется интерполяционной или аппроксимирующей.
Проще всего приближенное выражение производной получить с помощью соответствующей конечной разности.
Пользуясь определением производной функции, запишем:
f
′(x) = dy/dx ≈
h
x
f
h
x
f
)
(
)
(
−
+
,
(9.1) где h – конечная величина. Числитель представляет собой конечную разность первого порядка (1.1), следовательно,
f
′(x) ≈ ∆
1
y
/ h.

(9.2)
Таким образом, для получения приближенного значения производной в точке x необходимо вычислить два значения исходной функции f(x) – в точках x и x + h. Аналогично можно приближенное значение производной f
′(x) в точке x выразить с помощью значений исходной функции f(x) в точках x– h и x:
f
′(x) ≈
h
h
x
f
x
f
)
(
)
(
−
−

(9.3)
Представление производной f
′(x) формулами (9.1) или (9.3) означает, что на интервале [x, x + h] или [x – h, x] исходная функция f (x) заменяется линейной.
Более точные численные оценки производной функции можно получить, используя замену исходной функции f(x) интерполяционными полиномами.

135
Рассмотрим ситуацию, когда в нашем распоряжении имеется конечное число значений исследуемой функции f(x) в равноотстоящих узлах:
x
i
= a + i h, y
i
= f (x
i
) , i = 0, 1, …, n.
(9.4)
Теперь заменим исходную функцию f(x) первым полиномом
Ньютона (2.21). Для дальнейших преобразований полезно раскрыть скобки в числителях всех слагаемых:
P
n
(x) = y
0
+ q
∆
1
y
0
+
2 2
q
q
− ⋅
∆
2
y
0
+
6 2
3 2
3
q
q
q
+
−
∆
3
y
0
+
+
24 6
11 6
2 3
4
q
q
q
q
−
+
−
∆
4
y
0
+
120 24 50 35 10 2
3 4
5
q
q
q
q
q
+
−
+
−
∆
5
y
0
+…, (9.5) где q = (x
−
x
0
) / h.
Производную функции y = f(x) по аргументу x представим как производную сложной функции:
dx
dy
=
dq
dy
dx
dq
=
h
1
dq
dy

(9.6)
Дифференцирование полинома (9.5) по переменной x, с учетом (9.6), дает формулу для приближенного значения производной:
f
′(x) ≈ [∆
1
y
0
+
2 1
2
−
q
∆
2
y
0
+
6 2
6 3
2
+
−
q
q
∆
3
y
0
+
12 3
11 9
2 2
3
−
+
−
q
q
q
∆
4
y
0
+
+
120 24 100 105 40 5
2 3
4
+
−
+
−
q
q
q
q
∆
5
y
0
+ …] / h.
(9.7)
Если продифференцировать по переменной x полином (9.7), то мы получим формулу для приближенного значения второй производной:

136
f
′′(x) ≈ [ ∆
2
y
0
+ (q – 1)
∆
3
y
0
+
12 11 18 6
2
+
−
q
q
∆
4
y
0
+
+
12 10 21 12 2
2 3
−
+
−
q
q
q
∆
5
y
0
+ … ] / h
2

(9.8)
Напомним, что все конечные разности в формулах (9.2), (9.7) и (9.8) должны вычисляться для аргумента x = x
0
. При практическом использовании формул (9.5) и (9.6) в качестве узла x
0
выбирается узел таблицы, ближайший слева к значению аргумента x.
Если требуется вычислить значение первой или второй производной в каком-либо узле таблицы (9.2), то для этого достаточно в выражениях
(9.5) и (9.6) положить q = 0. Тогда получим следующие полезные формулы численного дифференцирования:
f
′(x
0
)
≈ [ ∆
1
y
0
− ∆
2
y
0
/ 2 +
∆
3
y
0
/ 3
− ∆
4
y
0
/ 4 +
∆
5
y
0
/ 5 …] / h.
(9.9)
f
′′(x
0
)
≈ [ ∆
2
y
0
–
∆
3
y
0
+
12 11
∆
4
y
0
–
6 5
∆
5
y
0
+…] / h
2

(9.10)
В выражениях (9.5) и (9.7) – (9.10) приведены несколько первых членов полинома Ньютона и его производных. В зависимости от числа узлов таблицы (9.5) и требуемой точности количество слагаемых в (9.7)
– (9.10) может быть уменьшено или увеличено. Заметим, что конечные разности
∆
i
y
0
в формулах (9.9) и (9.10) должны вычисляться для точки
x
0
, в которой ищутся производные.
Видно, что оценка производной (9.2) является частным случаем формулы (9.7), если степень полинома Ньютона положить равной единице (т.е. заменить неизвестную функцию f(x) линейной).
С помощью второго полинома Ньютона (2.26) возможно построить формулы численного дифференцирования, аналогичные вышеприведенным (9.7) – (9.10). Также в качестве дифференцируемых функций можно брать интерполяционные полиномы Стирлинга (2.31) и
Бесселя (2.32).
Погрешности формул (9.5) – (9.8) определяются абсолютной величиной первого из отбрасываемых членов. Так как конечная разность n-го порядка
∆
n
y
h
n
, то ясно, что погрешности оценок производных (9.1) и (9.3) имеют порядок
h, т.е. являются весьма грубыми приближениями.

137
Пример 1. Некоторая функция f (x) задана таблицей значений 9.1.
Таблица 9.1
i
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7
i
0 0,4 992 0,9 933 1,4 776 1,9 471 2,3 971 2,8 232 3,2 211
Пусть требуется вычислить первую и вторую производную данной функции f (x) в точке x = 0,24. Видно, что функция f (x) задана таблицей с равноотстоящими узлами, шаг таблицы h = 0,1. Узел 0,2 является ближайшим слева к заданному значению аргумента x = 0,24 и поэтому выбирается за x
0
. Для вычисления производных используем формулы
(9.7) и (9.8), построенные на основе первого интерполяционного полинома Ньютона.
Предварительно рассчитаем необходимые конечные разности пяти порядков в точке x = x
0
и сведем их в таблицу 9.2.
Таблица 9.2
x
i

y
i

∆y
i
∆
2
y
i
∆
3
y
i
∆
4
y
i
∆
5
y
i
0
,2 0,99 33 0,48 43
-
0,0148
-
0,0047 0,00 019 4,5
⋅
10
-5 0
,3 1,47 76 0,46 95
-
0,0195
-
0,0045 0,00 024 0
,4 1,94 71 0,45 00
-
0,0240
-
0,0043

0
,5 2,39 71 0,42 61
-
0,0282

0
,6 2,82 32 0,39 79


0
,7 3,22 11



Для заданного аргумента x = 0,24 величина q = (x
−
x
0
) / h = 0,4.
Расчет по формуле (9.7) с использованием только двух слагаемых дает величину 4,85726. Погрешность полученного значения определяется

138 абсолютной величиной третьего слагаемого формулы (9.7), которая равна 6,25 10
−4
. Следовательно, производная f
′(x = 0,24) равна 4,857 с тремя верными знаками после запятой.
Для получения более точного значения следует в формуле (9.7) использовать большее количество слагаемых. Например, расчет четырех слагаемых дает значение f
′(x = 0,24) = 4,8567 с четырьмя верными знаками после запятой. Учет следующих слагаемых в (9.7) не даст достоверных результатов, хотя бы потому, что исходные данные содержат ограниченное количество верных значащих цифр.
Расчет значения второй производной в точке x = 0,24 проводится по формуле (9.8) с использованием конечных разностей той же таблицы
9.2. Вычисление четырех слагаемых дает значение f
′′(x = 0,24) = –1,188 с тремя верными знаками после запятой.
Формулы, полученные с помощью интерполяционных полиномов
Ньютона, обладают существенным недостатком. Значение производной в определенной точке x
0
вычисляется по формулам (9.7) – (9.10) с использованием значений исходной функции f (x) только для аргументов x
≥x
0
. Иначе говоря, оценка численного значения производной рассчитывается на основе информации об исходной функции f (x) в области аргументов, превышающих значение x
0
Аналогично, формулы, полученные из второго интерполяционного полинома Ньютона, также страдают односторонностью оценки, так как используют значения исходной функции f (x) в точках, расположенных левее заданного аргумента. Большую точность дают более
«симметричные» формулы численного дифференцирования, которые используют значения заданной функции f (x) в точках x > x
0
и в x < x
0
, где x
0
– аргумент, для которого требуется отыскать значение производной. Такие формулы рассматриваются в следующем параграфе.

§ 9.2. Использование разложения в ряд Тейлора
Для практических целей наиболее целесообразным является представление производных формулами, использующими значения дифференцируемой функции в точках, расположенных симметрично относительно точки, в которой должна быть вычислена производная.
Такие формулы характеризуются высокой точностью, что позволяет

139 использовать сравнительно небольшое количество узлов для вычисления значений исходной функции f (x).
Получить
«симметричные» формулы численного дифференцирования можно, используя полиномы Стирлинга или
Лагранжа для набора равноотстоящих узлов. Более удобным является использование разложения данной функции f (x) в ряд Тейлора.
Пусть требуется найти значение производной заданной функции f (x) в точке x = x
0
. Будем полагать, что в нашем распоряжении имеется алгоритм, позволяющий получать значения функции f (x) в некотором конечном интервале аргументов x в окрестности точки x
0
Вычислим значения данной функции f (x) в равноотстоящих точках
x
i
= x
0
+ i h,
i
= 0,
±1, ±2, ±3, … .
(9.11)
Величину h выберем малой. Тогда точки x
1
= x
0
+ h и x
−1
= x
0
− h можно полагать находящимися в окрестности точки x
0
. Запишем следующие разложения в ряд Тейлора для функции f (x):
f
(x
1
) = f (x
0
) + h f
′(x
0
) +
2 2
h
f
′′(x
0
) +
!
3 3
h
f
′′′(x
0
) + …,
(9.12)
f
(x
−1
) = f (x
0
)
− h f ′(x
0
) +
2 2
h
f
′′(x
0
)
−
!
3 3
h
f
′′′(x
0
) + … .
(9.13)
Сначала в разложениях (9.12) и (9.13) отбросим все члены 3-го порядка и выше. Из-за малости величины h такое приближение вполне допустимо. После этого вычтем выражение (9.13) из (9.12) и получим следующее уравнение:
f
(x
1
)
− f (x
−1
) = 2 h f
′(x
0
).
Отсюда сразу получается формула для приближенного вычисления производной в точке x = x
0
:
f
′(x
0
)
≈
h
h
x
f
h
x
f
2
)
(
)
(
0 0
−
−
+

(9.14)

140
Погрешность последней формулы определяется величиной первого отброшенного члена разложений (9.12) и (9.13) и оценивается величиной h
2
f
′′′(x
0
) / 6 , т.е. имеет порядок
h
2
. Для сравнения отметим, что погрешность формулы (9.1) на порядок больше.
Формула (9.14) называется 3-точечной оценкой первой производной функции f (x), хотя для получения величины f
′(x
0
) достаточно вычислить значения исходной функции f (x) только в двух точках, симметричных относительно аргумента x
0
Оценка производной (9.14) может быть уточнена. Для этого в разложениях (9.12) и (9.13) отбросим все члены, начиная с 4-го. После этого проведем вычитание (9.13) из (9.12) и перегруппируем члены. В результате получится уравнение, связывающее значения функции f (x) и ее производных в разных точках:
f
(x
1
)
− f (x
−1
) = 2 h f
′(x
0
) + h
3
f
′′′(x
0
) / 3 .
(9.15)
Теперь запишем разложения в ряд Тейлора для функции f (x) в точках x
2
= x
0
+ 2 h и x
−2
= x
0
− 2 h, ограничиваясь членами 3-го порядка включительно:
f
(x
2
) = f (x
0
) + 2 h f
′(x
0
) + 4 2
2
h
f
′′(x
0
) + 8
!
3 3
h
f
′′′(x
0
) + …,
(9.16)
f
(x
−2
) = f (x
0
)
− 2 h f ′(x
0
) + 4 2
2
h
f
′′(x
0
)
− 8
!
3 3
h
f
′′′(x
0
) + ….
(9.17)
Проведем вычитание разложения (9.17) из (9.16):
f
(x
2
)
− f (x
−2
) = 4 h f
′(x
0
) + 8 h
3
f
′′′(x
0
) / 3.
Из последнего уравнения выразим h
3
f
′′′(x
0
) и подставим в (9.15).
При этом получится уравнение, которое можно разрешить относительно величины f
′(x
0
). Таким образом, получается еще одна приближенная формула для вычисления первой производной функции f (x) в точке
x
=x
0
:

141
f
′(x
0
)
≈
h
h
x
f
h
x
f
h
x
f
h
x
f
12
)
2
(
)
(
8
)
(
8
)
2
(
0 0
0 0
+
−
+
+
−
−
−
(9.18)
Погрешность формулы (9.18) оценивается величиной h
5
f
(V)
(x
0
) / 120.
Следовательно, погрешность формулы (9.18) на 3 порядка меньше погрешности оценки (9.14).
Формула (9.18) называется 5-точечной оценкой первой производной функции f (x). Заметим, что вычисления значения функции f (x) в точке
x
=x
0
не требуется.
Перейдем к выводу формул численного дифференцирования для вторых производных данной функции. Для этого вернемся к разложениям (9.12) и (9.13), отбросим все члены выше 3-го порядка и сложим их. Получим уравнение:
f
(x
1
) + f (x
−1
) = 2 f (x
0
) + h
2
f
′′(x
0
), из которого непосредственно получается формула приближенного вычисления второй производной в точке x = x
0
:
f
′′(x
0
)
≈
2 0
0 0
)
(
)
(
2
)
(
h
h
x
f
x
f
h
x
f
−
+
−
+

(9.19)
Эта формула справедливо называется 3-точечной оценкой второй производной, так как величина f
′′(x
0
) выражается через значения исходной функции f (x) в трех точках, симметричных относительно узла
x
0
Погрешность оценки (9.19) также определяется отбрасываемыми членами разложения и приблизительно равна по абсолютной величине
h
2

| f
(IV)
(x
0
)
| / 12, т.е. имеет порядок
h
2
Для нахождения более точной оценки второй производной необходимо воспользоваться разложениями в ряд Тейлора исходной функции до 4-го члена включительно в точках x
1
= x
0
+ h, x
−1
= x
0
− h, x
2
= x
0
+ 2h и x
−2
= x
0
− 2h.
f
(x
1
) = f(x
0
) + h f
′(x
0
) +
2 2
h
f
′′(x
0
) +
!
3 3
h
f
′′′(x
0
) +
!
4 4
h
f
(IV)
(x
0
),
(9.20)

142
f
(x
−1
) = f(x
0
)
− h f ′(x
0
) +
2 2
h
f
′′(x
0
)
−
!
3 3
h
f
′′′(x
0
) +
!
4 4
h
f
(IV)
(x
0
),
(9.21)
f
(x
2
) = f(x
0
) + 2h f
′(x
0
) +
2 4
2
h
f
′′(x
0
) +
!
3 8
3
h
f
′′′(x
0
) +
!
4 16 4
h
f
(IV)
(x
0
),
(9.22)
f
(x
−2
) = f(x
0
)
− 2h f ′(x
0
) +
2 4
2
h
f
′′(x
0
)
−
!
3 8
3
h
f
′′′(x
0
) +
!
4 16 4
h
f
(IV)
(x
0
).
(9.23)
Сложим попарно разложения (9.20) и (9.21), а также (9.22) и (9.23).
Получим уравнения:
f
(x
1
) + f (x
−1
) = 2 f (x
0
) + h
2
f
′′(x
0
) +
12 4
h
f
(IV)
(x
0
),
(9.24)
f
(x
2
) + f (x
−2
) = 2 f (x
0
) + 4 h
2
f
′′(x
0
) +
3 4
4
h
f
(IV)
(x
0
).
(9.25)
Выразим произведение h
4
f
(IV)
(x
0
) из (9.24) и подставим (9.25). Таким путем получим соотношение, связывающее значения функции f (x) в пяти разных узлах и второй производной в точке x = x
0
:
f
(x
2
) + f (x
−2
) = 16 [f (x
1
) + f (x
−1
)]
− 30 f (x
0
)
− 12 h
2
f
′′(x
0
).
Отсюда получается новая формула для приближенного вычисления второй производной заданной функции в точке x = x
0
:
f
′′(x
0
)
≈
≈
2 0
0 0
0 0
12
)
2
(
)
(
16
)
(
30
)
(
16
)
2
(
h
h
x
f
h
x
f
x
f
h
x
f
h
x
f
+
−
+
+
−
−
+
−
−
(9.26)

143
Формула (9.26) называется 5-точечной оценкой второй производной, так как использует пять значений исходной функции f (x) в пяти равноотстоящих точках. Численное дифференцирование по формуле дает погрешность (9.26) приблизительно равную h
4

| f
(VI)
(x
0
)
| / 90, что на два порядка меньше погрешности оценки (9.19).
Пример 2. Вновь рассмотрим функцию, заданную таблицей 9.1 в примере 1 этого параграфа. Вычислим первую и вторую производные в узле таблицы x = 0,3, приняв это значение за центральный узел x
0
Используя величину шага h = 0,1, получим по формуле (9.14) значение первой производной f
′(x = 0,3) = 4,769, а по формуле (9.18) − значение
f
′(x = 0,3) = 4,777. Расчеты второй производной по формулам (9.19) и
(9.26) дадут значения
−1,48 и −1,4825 соответственно.
Для сравнения вычислим обе производные в той же точке x = 0,3 с помощью формул (9.9) и (9.10), полученных из первого интерполяционного полинома Ньютона. Напомним, что при этом должны использоваться конечные разности, рассчитанные в точке x =
0,3, т.е. расположенные во второй строке данных таблицы 9.2.
Использование двух слагаемых в сумме (9.9) дает результат f
′(x = 0,3) =
4,792, а четырех слагаемых – значение f
′(x = 0,3) = 4,777. Расчеты второй производной по формуле (9.10) с использованием двух и четырех слагаемых дают значения
−1,5 и −1,474 соответственно.
Теперь можно открыть маленький секрет. В таблице 9.1 расположены значения функции y = 5sin(x). Аналитические расчеты первой и второй производных для аргумента x = 0,3 дают значения f
′(x)
= 4,7767 и f
′′(x) = −1,4776. Сравнение результатов, полученных разными методами, позволяет сделать следующие выводы. Расчет величины f
′(x) по 5- точечной формуле (9.18) и с помощью выражения
(9.9) с четырьмя слагаемыми дают 3 верные цифры после запятой, что для многих практических задач вполне достаточно. Расчеты по 3- точечной формуле (9.14) и по формуле (9.9) с двумя слагаемыми дают менее точные результаты.
Вычисления второй производной f
′′(x) по 5-точечной формуле (9.26) и по формуле (9.10) с использованием четырех слагаемых дают результат с точностью до сотых долей.
Таким образом, оба подхода к вычислению производных, изложенные в § 9.1, 9.2, позволяют достичь определенной точности, однако расчеты по формулам § 9.2 значительно короче. С другой стороны, формулы § 9.1 необходимы для случаев, когда требуется вычислять производные в точках, не совпадающих с узлами.

144
Обращает на себя внимание более низкая точность вычислений второй производной (по сравнению с первой).
Естественным способом понижения погрешности численных расчетов величин производных кажется уменьшение шага h таблицы узлов (в тех случаях, когда это возможно). Однако принципиальная некорректность процедуры численного дифференцирования, о которой упоминалось в начале главы, приводит к иному результату. Проверка показывает [9], что с уменьшением шага h таблицы узлов точность численного дифференцирования сначала возрастает, затем начинает существенно убывать.
Обратим внимание на то, что в знаменателях всех формул численного дифференцирования стоит малое число, а это приводит к погрешностям округления при компьютерных вычислениях. При заданном шаге h относительная погрешность вычисления для второй производной больше, чем для первой.
На практике при численном расчете производных следует ограничиваться минимальным количеством значащих цифр.

§ 9.3. Численное дифференцирование при произвольном расположении узлов
Иногда требуется получить оценку значения производной k-го порядка, когда вся информация о дифференцируемой функции сводится к таблице ее значений y
i
в неравноотстоящих узлах x
i

(i = 0, 1, ... , n). В этих случаях по заданной таблице можно построить полином Лагранжа, а затем его аналитически дифференцировать. Однако это сопряжено с достаточно громоздкими алгебраическими преобразованиями.
Альтернативой в такой ситуации является численное дифференцирование с использованием метода

Скачать 0.96 Mb.

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