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

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(x0 +qh) = y(x0) + qDy0 + + + + …+ , h = x i+1 - x i; i = 1, 2, …, n.

(4)

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

y(x) » y(x₀)+ qDy₀ + + + + … (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(x0 + 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. Отметим, что в этом случае для вычисления значений конечных разностей Dy_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₀)+ qDy₀, проходящей через точки x₀ и x₁, а её производная (6) приближённо заменяется выражением: , где Dy₀ = y(x₁) – y(x₀)

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

(17)

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

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

P₂(x₀ +qh) = y(x₀)+ qDy₀ + ,

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

, где Dy₀ = y(x₁) – y(x₀); D²y₀ = y(x₂) – 2y(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₀) + qDy₀ + + ,

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

,

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

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

(21)

(22)

(23)

(24)