Конечноразностные формулы численного дифференцирования на основе интерполяционных полиномов

2.1 Построение формул численного дифференцирования на основе первой интерполяционной формулы Ньютона.

Пусть дифференцируемая функция y = f (x) задана своими значениями
y ₀, y ₁, y ₂, …, y_n, на сетке равноотстоящих узлов : x ₀, x ₁, x ₂, …, x_n отрезка [ a, b ], т.е. в данном случае:

x_i = x ₀ + ih, где i = 1, 2, …, n,

где h > 0 - шаг сетки (или таблицы).

Для нахождения на отрезке [ a, b ] значений производных y ¢ = f¢ (x), y ¢¢ = f¢¢ (x) и т.д., функцию y = f (x) в узлах отрезка [ a, b ] (приближённо) заменяют первым интерполяционным полиномом Ньютона, построенным для данной системы равноотстоящих узлов .

y (x)» Pn (x 0 + qh) = y (x 0) + q D y 0 + + + + …+ , h = x i+ 1 - x i; i = 1, 2, …, n.

(4)

Производя перемножение биномов в формуле (4) получим:

y (x)» y (x ₀) + q D y ₀ + + + + … (5)

Поскольку , то в соответствии с выражениями (4) и (5) функцию P_n (x ₀ + qh) можно рассматривать как сложную функцию переменной x. Поэтому в соответствии с правилом дифференцирования сложной функции имеем:

f¢_x (q (x)) @ [ P_n (x ₀ + q (x) h)] ¢_x = [ P_n (x ₀ + q (x) h)] ¢_q × q¢_x, где q¢_x

и, выполняя дифференцирование формулы (5), получаем:

yx ¢(x) = y¢q (x 0 + q (x) h)× q¢x @

(6)

Далее по аналогии с выражением (6) для y ¢(x) можно написать выражение для y ²(x) и y ¢¢¢(x):

y ¢¢(x) @ (7)

y ¢¢¢(x) @ (8)

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

Наиболее важной в приложениях является простейшая формула для аппроксимация второй производной с помощью разностного отношения:

y ¢¢(x)» . (9)

Формула (9) получается из выражения (5) при n = 2, на промежутке x Î(x ₀, x ₂).

Аппроксимация второй производной в виде разностного отношения (9) часто востребована при построении конечноразностных методов решения краевых задач для обыкновенных дифференциальных уравнений второго порядка и при разработке численных методов решения уравнений в частных производных.

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

2.2 Численное дифференцирование в узловых точках таблицы.

Иногда требуется находить производные функции y (x) в самих табличных точках x _i (т.е. узлах сетки). В этом случае формулы численного дифференцирования упрощаются.

При построении интерполяционного полинома Ньютона каждое табличное значение x_i, в котором нужно вычислить значение производной f¢ (x_i), может быть принято за базовый узел x ₀. Поэтому в случае, когда базовым узлом x ₀ является узел x_i, переменная в формулах (6) – (8) будет равна и формулы (6) – (8) перепишутся в виде:

y ¢(x) @ ; y ¢¢(x) @ ; y ¢¢¢(x) @ ; где .

(10)

Формулами (10) задаются выражения для аппроксимации производных y ¢(x), y ¢¢(x), y ¢¢¢(x) функции y = f (x) вблизи табличной точки x _i, y (x_i).

При подстановке в формулы (10) вместо x значения базового узла x_i
(x = x_i) переменная обращается в нуль и для значений производных в узлах таблицы, т.е. в точках x = x _i получим следующие выражения:

y ¢(x_i) @ (11)

y ¢¢(x_i) @ (12)

y ¢¢¢(x_i) @ (13)

Чтобы построить аппроксимации (11) - (13) для y ¢(x_i) и y ¢¢(x_i), включающие, например, конечные разности пятого порядка, нужно по формулам (4) или (5) построить первый интерполяционный полином Ньютона пятого порядка (n = 5), выбрав в качестве базового узла точку x_i. Отметим, что в этом случае для вычисления значений конечных разностей D y_i, D² y_i, D³ y_i, D⁴ y_i, D⁵ y_i необходимо использовать шесть строк из таблицы значений функции y = f (x):

y (x_i), ® y (x_{i+ 1}), ® y (x_{i+ 2}), ® y (x_{i+ 3}), ® y (x_{i+ 4}), ® y (x_{i+ 5}).

Формулы численного дифференцирования (10) - (13), построенные на основе первой интерполяционной формулы Ньютона, более эффективны для использования в начальной части таблицы значений функции, т.е. когда в таблице значений функции «ниже» базового узла x_i имеется достаточное количество строк, необходимых для вычисления соответствующих значений конечных разностей.

Аналогичные формулы численного дифференцирования могут быть получены и на основе второй интерполяционной формулы Ньютона. Они более эффективны при выборе базовых узлов в конечной части таблицы значений функции. Также можно вывести и ряд «центральных» формул численного дифференцирования на основе рассмотренных ранее интерполяционных формул Гаусса, Бесселя, Стирлинга, такие формулы более эффективно «работают» в близи центральных узловых точек таблицы значений функции [1, стр. 567 - 570].

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

1. Для нахождения значений производных y ¢(x), y ¢¢(x), y ¢¢¢(x) … функции y (x) в некоторой фиксированной (межузловой) точке x^* таблицы значений функции
y = f (x) в качестве базового узла x_i в формулах (10) необходимо выбрать ближайшее к x^* табличное значение аргумента x_i (лежащее левее точки x^*).

2. При желании получить аппроксимации производных y ¢(x), y ¢¢(x), y ¢¢¢(x) … с наибольшей точностью максимальный порядок конечных разностей в формулах (10) - (13) определяется при построении таблицы конечных разностей функции y = f (x). При этом, исходя из анализа таблицы конечных разностей, необходимо оценить оптимальную степень n интерполирующего полинома P_n (x) и сохранить в формулах (10) - (13) конечные разности до D ⁿy_i включительно.

3. Отметим, что занижение степени интерполирующего полинома P_n (x) (4) - (5) отрицательно сказывается на точности полученных формул численного дифференцирования так же, как и необоснованное завышение его степени. Завышение степени полинома выше оптимальной приводит к учёту в формулах численного дифференцирования ошибок округления. Это происходит потому, что из-за ограниченной точности задания значений сеточной функции y (x) конечные разности высоких порядков начинают возрастать с увеличением их порядка вследствие «расползания» исходных ошибок округления и, следовательно, конечны разности высоких порядков могут содержать большие ошибки.

2.3 Связь приближённых значений производной функции с узловыми значениями самой функции.

Адаптируем полученную выше формулу численного дифференцирования:

y ¢(x) @ (6)

для наиболее часто используемых практических ситуаций.

Если показатель - n степени интерполяционного полинома Ньютона (5), лежащего в основе формулы (6), зафиксировать соответственно равным: n = 1, n = 2, n = 3, то это будет соответствовать ситуации, когда дифференцируемая функция y = f (x) на отрезке интерполирования заменяется соответственно:

a. прямой (n = 1) - случай линейной интерполяции по двум точкам x ₀; x ₁;

b. параболой (n = 2) - случай квадратичной интерполяции по трём точкам x ₀; x ₁; x ₂;

c. кубической параболой (n = 3) - случай кубической интерполяции по четырём точкам x ₀; x ₁; x ₂; x ₃.

Для рассматриваемых случаев из формулы (6) получаются следующие выражения для приближённых значений производной y ¢(x) функции y = f (x):

для x Î[ x ₀; x ₁] (14)

– оценка (аппроксимация) производной на основе линейной интерполяции.

для x Î[ x ₀; x ₂] (15)

– оценка (аппроксимация) производной на основе квадратичной интерполяции.

для x Î[ x ₀; x ₃] (16)

– оценка (аппроксимация) производной на основе кубической интерполяции.

Получим на основе формул (14) – (16) аппроксимации для y ¢(x) = f¢ (x) в узлах интерполяции, выраженные через значения функции y = f (x) в узлах сетки. Для этого учтём, что при подстановке в формулы (14) - (16) вместо переменной x узловых значений x ₀, x ₁, x ₂; x ₃ переменная будет соответственно равна: 0, 1, 2, 3.

a) Случай линейной интерполяции: n = 1, x Î[ x ₀; x ₁].

В этом случае на отрезке [ x ₀; x ₁] функция y = f (x) заменяется прямой линией P ₁(x ₀ + qh) = y (x ₀) + q D y ₀, проходящей через точки x ₀ и x ₁, а её производная (6) приближённо заменяется выражением: , где D y ₀ = y (x ₁) – y (x ₀)

Поэтому из (6) в данном случае имеем:

(17)

б) Случай квадратичной интерполяции: n = 2, x Î[ x ₀; x ₂].

В этом случае на отрезке [ x ₀; x ₂] функция y = f (x) заменяется параболой:

P ₂(x ₀ + qh) = y (x ₀) + q D y ₀ + ,

проходящей через точки x ₀, x ₁, x ₂, а её производная f¢ (x) (6) заменяется выражением:

, где D y ₀ = y (x ₁) – y (x ₀); D² y ₀ = y (x ₂) – 2 y (x ₁) + y (x ₀);

Откуда, учитывая, что при x = x ₀ q = 0, а при x = x ₁ q = 1, при x = x ₂ q = 2, получаем:

(18)

(19)

(20)

в) Случай кубической интерполяции: n = 3, x Î[ x ₀; x ₃].

В этом случае на отрезке [ x ₀; x ₃] функция y = f (x) заменяется кубической параболой, проходящей через точки x ₀, x ₁, x ₂, x ₃:

P ₃(x ₀ + qh) = y (x ₀) + q D y ₀ + + ,

а её f¢ (x) (6) заменяется выражением:

,

где D y ₀ = y (x ₁) – y (x ₀); D² y ₀ = y (x ₂) – 2 y (x ₁) + y (x ₀); D³ y ₀ = y (x ₃) – 3 y (x ₂) +3 y (x ₁) – y (x ₀);

Так как при x = x ₀ q = 0, при x = x ₁ q = 1, при x = x ₂ q = 2, при x = x ₃ q = 3, получаем:

(21)

(22)

(23)

(24)