Численное дифференцирование

Пусть на отрезке [a, b] в узлах одномерной сетки

h_x ={x_i / x_i = x_i–1 + h_i, h_i > 0, i = 1, 2, 3, …, n; x₀ = a, x_n = b}

заданы y_i = f(x_i), i = 0, 1, 2, …, n – значения дифференцируемой на этом отрезке функции f(x).

Для упрощения последующих записей будем считать, что сетка h_x равномерная.

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

(4.2.1)

Производные в точках a и b (граничных точках отрезка) могут быть в этом случае определены как односторонние пределы:

(4.2.2)

Здесь Dx ® +0 означает, что Dx стремится к нулю, оставаясь при этом больше нуля.

На основании (4.2.1) и (4.2.2) приближенные значения производных
в узлах сетки h_x могут быть вычислены с помощью левостороннего, правостороннего и центрального конечно-разностных отношений соответственно:

(4.2.3)

(4.2.4)

(4.2.5)

По формулам (4.2.3) и (4.2.4) определяют так называемые левостороннюю и правостороннюю разностные производные, по формуле (4.2.5) – центральную разностную производную.

Аналогично могут быть вычислены приближенные значения производных второго и более высоких порядков. Так, для определения производной второго порядка для i = 1, 2, …, n – 1 может быть использована следующая приближенная формула:

(4.2.6)

Пусть теперь d – абсолютная погрешность приведенных выше простейших формул численного дифференцирования. Здесь имеется в виду погрешность, возникающая при аппроксимации производной соответствующим конечно-разностным отношением, другими словами, погрешность метода. Для простоты при оценке погрешности аппроксимации сетку h_x будем считать равномерной, с шагом h. Тогда для (4.2.3) и (4.2.4) имеет место оценка

(4.2.7)

где M₂ – максимум модуля производной второго порядка функции f(x) на отрезке [a, b].

Согласно (4.2.7) формулы (4.2.3) и (4.2.4) имеют погрешность первого порядка малости по h и точны (d = 0) для полиномов до первого порядка включительно.

Для формулы (4.2.5) оценка погрешности

(4.2.8)

где M₃ – максимум модуля производной третьего порядка функции f(x) на отрезке [a, b].

Согласно (4.2.8) формула (4.2.5) имеет погрешность второго порядка малости по h и точна (d = 0) для полиномов до второго порядка включительно.

Из приведенных выше теоретических оценок погрешности метода для формул (4.2.3)–(4.2.5) следует, что при стремлении шага сетки h к нулю погрешность получаемого приближенного значения производной также будет стремиться к нулю. Однако на практике этого не происходит, напротив, при h ® 0 погрешность получаемых результатов численного дифференцирования может неограниченно возрастать. Это происходит потому, что в соотношениях (4.2.7) и (4.2.8) оценивается только та часть погреш-ности, которая имеет место из-за замены производной конечно-разностным отношением и называется погрешностью аппроксимации. Другая часть результирующей погрешности обусловлена погрешностью производимых вычислительных операций над имеющими ненулевую погрешность данными значениями функции f(x) в узлах сетки h_x.

Действительно, пусть абсолютная погрешность y_i, i = 0, 1, 2, …, n может быть оценена сверху числом s > 0. Тогда в результате вычисления конечно-разностного отношения, например, по формуле (4.2.3) будет
получен результат, имеющий абсолютную вычислительную погрешность, которая может быть оценена сверху числом 2s / h. Значение s не зависит от h, и при h ® 0 оценка 2s / h неограниченно возрастает. Аналогично может быть проанализировано поведение вычислительной погрешности при определении конечно-разностных отношений по формулам (4.2.4) и (4.2.5).

Другими словами, при наличии ненулевой погрешности данных y_i,
i = 0, 1, 2, …, n, представляющих собой значения функции f(x) в узлах сетки h_x, суммарная погрешность приближенного расчета производных (складывается из погрешности аппроксимации и вычислительной погрешности) с помощью конечно-разностных отношений (4.2.3) – (4.2.5) может неограниченно возрастать при стремлении шага сетки к нулю. В связи с отмеченными свойствами погрешности операция численного дифференцирования в вычислительном отношении является некорректной.

С практической точки зрения при вычислениях по формулам
(4.2.3) – (4.2.5) необходимо стремиться к тому, чтобы вычислительная погрешность не превосходила погрешности аппроксимации. Применительно к формуле (4.2.3) это означает, что должно выполняться неравенство

(2s / h) £ (M₂h / 2). (4.2.9)

Это неравенство может служить основанием для задания s – требуемой точности вычисления значений функции f(x) в узлах сетки h_x, если шаг h сетки задан, или для определения минимально допустимого значения h, если значение s не может быть изменено. Например, из (4.2.9) следует, что s не должно превышать значения M₂h² / 4.

В целом организация вычислений приближенного значения производной с использованием одной из формул (4.2.3) – (4.2.5) очевидна и не вызывает затруднений. Вместе с тем необходимо отметить, что эти простейшие формулы численного дифференцирования позволяют вычислять приближенное значение производной f(x) только в узлах сетки h_x и не позволяют производить вычисления в произвольной точке x^* Î (a, b), не совпадающей ни с одним из этих узлов. Единственное, что можно предпринять для произвольной точки x^* Î (a, b), это проинтерполировать по двум (или более) соседним с x^* узлам сетки h_x те приближенные значения производных, которые будут найдены в этих узлах по одной из формул (4.2.3) – (4.2.5). При этом нет смысла добиваться большей точности интерполирования, чем точность приближенного дифференцирования.

Имеется более общий подход к приближенному вычислению производных для таблично заданных функций, который заключается в использовании интерполяционного полинома. Действительно, если для заданной
в узлах сетки h_x функции f(x) известен интерполяционный полином

P_n(x) = a₀ + a₁x + a₂x² + … + a_nxⁿ, (4.2.10)

то, полагая, что функция f(x) m раз (m £ n) дифференцируема на отрезке [a, b], производную f^(m)(x^*) в любой точке x^* Î [a, b] можно приближенно вычислить как производную соответствующего порядка интерполяционного полинома P_n(x):

f^(m)(x^*) » P_n^(m)(x^*). (4.2.11)

Погрешность нахождения производной по формуле (4.2.11) определяется модулем производной соответствующего порядка от R_n(x^*) – остаточного члена интерполяции, задаваемого в соответствии с соотношением

f(x) = P_n(x) + R_n(x), x Î [a, b]. (4.2.12)

Важно отметить, что в общем случае из малости значений R_n(x)
не следует малость значений производных от R_n(x) и погрешность численного дифференцирования, как правило, превышает погрешность интерполяции.

Пусть на произвольном отрезке [x_{i – 1}, x_{i + 1}] Ì [a, b] по точкам x_{i – 1}, x_i
и x_{i + 1} построен интерполяционный полином в форме Лагранжа (формула (3.2.9)):

(4.2.13)

Для отыскания приближенного значения производной в произвольной точке x Î [x_{i – 1}, x_{i + 1}] в соответствии с (4.2.11) найдем выражение для первой производной P₂(x):

(4.2.14)

Если x = x_i, получим

(4.2.15)

В случае равномерной сетки формула (4.2.15) и формула центральной разностной производной (4.2.5) эквивалентны.

Теперь для приближенного вычисления производных воспользуемся интерполяционным полиномом в форме Ньютона (формула (3.3.1)). Обозначим через D^my_i табличные (конечные) разности m-го порядка в i-m узле сетки h_x:

(4.2.16)

Для упрощения записи у конечной разности первого порядка верхний индекс опускается.

Предположим, что сетка h_x равномерная, и введем обозначение

z: = (x – x₀) / h, x Î [a, b]. (4.2.17)

Тогда интерполяционный полином в форме Ньютона можно представить как

(4.2.18)

а приближенное равенство (4.2.11) для первой производной перепишется (для полинома, построенного по пяти узлам) в виде

(4.2.19)

Если x совпадает с одним из узлов сетки h_x, то этот узел всегда можно принять за x₀. Тогда x = x₀, z = 0 и формула (4.2.19) существенно упростится. Перепишем (4.2.19) для x = x₀ в общем случае при интерполировании по (n + 1) узлу сетки:

(4.2.20)

Очевидно, что при интерполировании по двум узлам (n = 1) фор-мула (4.2.20) эквивалентна формуле правосторонней разностной производной (4.2.4).

Погрешность аппроксимации производной функции f^(m)(x) производной интерполяционного полинома P_n^(m)(x) зависит от m-порядка производной и n-порядка полинома.

Можно показать, что при аппроксимации первой производной функции f(x) первой производной интерполяционного полинома, построенного по двум узлам, погрешность аппроксимации будет первого порядка по h. Если использовать интерполяционный полином, построенный по трем узлам, то погрешность аппроксимации будет второго порядка по h.

Помимо погрешности аппроксимации заметное влияние на погрешность получаемых приближенных значений производной оказывает вычислительная погрешность. В частности, известно, что погрешность вычисления конечных разностей быстро увеличивается с ростом их порядка. Поэтому при расчетах с применением формул (4.2.19) и (4.2.20) в общем случае не следует использовать большое количество слагаемых, целесообразно ограничиваться лишь двумя, например:

(4.2.21)

Для второй производной также можно использовать приближенное равенство

(4.2.22)

Пусть теперь заданы точка x^* Î (a, b) и значения функции f(x) на отрезке [a, b] в узлах одномерной сетки h_x. Для упрощения записи будем считать, что x^* не совпадает ни с одним из узлов сетки. На основании проведенного анализа процесс приближенного вычисления производной функции f(x) в точке x^* Î (a,b) с использованием интерполяционного полинома в форме Ньютона (оставляя два слагаемых в формуле (4.2.19)) можно представить, например, в виде следующего алгоритма:

1. На отрезке [a, b] задать одномерную равномерную сетку

h_x = {x_i / x_i = x_{i – 1} + h, h > 0, i = 1, 2, 3, …, n; x₀ = a, x_n = b}

и значения y_i = f(x_i) в узлах сетки x_i, i = 0, 1, 2, …, n.

Задать x^* Î (a, b), x^* ¹ x_i, i = 0, 1, 2, …, n.

2. Определить узел x_k сетки h_x исходя из условий x^* Î (x_k, x_{k + 1});
0 £ k < n – 1.

3. Если x_k найден, следует узлы x_k, x_{k + 1}, x_{k + 2} сетки h_x обозначить как x₀, x₁, x₂ соответственно и перейти к выполнению п. 5. Иначе – к выполнению п. 4.

4. Положить h: = –h. Узлы x_n, x_{n – 1}, x_{n – 2} сетки h_x обозначить как x₀, x₁, x₂ соответственно.

5. Вычислить z: = (x^* – x₀) / h.

6. Вычислить конечные разности Dy₀ и D²y₀.

7. Вычислить приближенное значение f¢(x) по формуле

8. Процесс завершен.

Действия, выполненные в соответствии с п. 1 приведенного алгоритма, позволяют определить узел сетки, ближайший к целевой точке x^* и находящийся на оси Ох левее этой точки. Выполнение условия 0 £ k < n – 1 обеспечивает возможность вычисления конечной разности второго порядка D²y₀. Если такой узел x_k не буден найден, это будет означать, что точка x^* находится в нижней части таблицы и для решения задачи необходимо воспользоваться интерполяционным полиномом в форме Ньютона для интерполирования назад (п. 3).

В заключение дадим ряд практических рекомендаций по выполнению численного дифференцирования.

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

Таблица 4.2.1

2. Процедура численного дифференцирования является с вычислительной точки зрения некорректной. Для обеспечения вычислительной
устойчивости необходимо ограничивать снизу величину шага h сетки h_x, не допуская ее неоправданного уменьшения.

3. На практике производные третьего и более высоких порядков численно определяются крайне редко из-за большой вычислительной погрешности получаемых результатов.

4. Использование интерполяционного полинома для реализации процедуры численного дифференцирования позволяет в общем случае достигать меньшей погрешности аппроксимации, чем применение простейших разностных отношений вида (4.2.3) – (4.2.5). Однако порядок интерполяционного полинома желательно ограничивать, не допуская применения полиномов более 5-го или 6-го порядков.