Интерполяция и сплайны

(главы из пособия, ч. 2 В.И.Середы)

 

Глава 3. Интерполирование и аппроксимация функций, заданных таблично

§ 3.1. Основные понятия

В данной главе рассматривается случай табличного задания функции одной переменной. С математической точки зрения это означает, что значения функции f(x) на промежутке [a, b] заданы в узлах одномерной сетки

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},         (3.1.1)

где xi, i = 0, 1, 2, …, n – узлы одномерной сетки hx; hi, i = 1, 2, 3, …, n – шаг сетки hx. Если hi = h " i = 1, 2, 3, …, n, сетка hx равномерная.

В случае табличного задания функции на отрезке [a, b] часто возникает проблема восстановления значений функции f(x) в произвольной точке x Î (a, b) между узлами сетки hx. Такую задачу называют задачей интерполирования (или интерполяции) значений функции в точку x Î (a, b).

Наряду с задачей интерполяции значений функции решается задача определения возможных значений таблично заданной на отрезке [a, b] функции f(x) в точках, не принадлежащих этому отрезку. Эту задачу также называют задачей экстраполяции значений функции в точку x Ï [a, b].

В более общей постановке ставится задача замены таблично заданной на отрезке [a, b] функции f(x) некоторой функцией j (x) из заданного класса функций (например, полиномов той или иной степени, экспоненциальных функций и т. д.). При этом выбор функции j (x) из заданного класса осуществляется таким образом, чтобы ее значения были в наибольшей степени согласованы в смысле некоторого критерия с табличными данными. Такая задача называется задачей аппроксимации или приближения таблично заданных функций. Формально задачи интерполяции и экстраполяции таблично заданной функции решаются посредством аппроксимации ее функциями заданного класса.

Рассмотрим некоторые методы решения задач интерполяции.

Приведем формальную постановку задачи полиномиальной интерполяции таблично заданной на отрезке [a, b] функции f(x). Пусть на отрезке [a, b] задана одномерная сетка

hx = {xi / xi = x 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).

Функция j (x) называется интерполирующей или интерполяционной для f(x) на [a, b], если ее значения j (x0), j (x1), …, j (xn) в узлах сетки
x0, x1, x2,…, xn совпадают с заданными значениями yi функции f(xi), i = 0, 1, 2, …, n. Геометрически это означает, что график функции j (x) проходит так, что по крайней мере в точке (n + 1) (в точках, соответствующих узлам сетки hx) на отрезке [a, b] он пересекает или касается истинного графика функции f(x).

Пусть в качестве интерполяционной функции j (x) выбран полином Pn(x) – алгебраический многочлен степени n:

Pn(x) = a0 + a1x + a2x2 + … + anxn.

Тогда задача полиномиальной интерполяции формулируется так: для функции f(x), заданной таблично, определить полином Pn(x), при котором выполняются условия

Pn(xi) = yi , i = 0, 1, 2, …, n.                                                                (3.1.2)

Для того чтобы определить полином, необходимо найти значения
его (
n + 1) коэффициента: a0, a1, a2, …, an. Для этого имеется (n + 1) условие (3.1.2). Иначе говоря, для того чтобы полином был интерполяционным для таблично заданной функции f(x), достаточно, чтобы его коэффициенты удовлетворяли системе уравнений

                                                           (3.1.3)

Поскольку все узлы сетки hx различны (hi > 0 "i), то определитель этой системы (определитель Вандермонда) отличен от нуля, следовательно, решение системы (3.1.3) существует и является единственным. Таким образом, чтобы построить интерполяционный полином Pn(x), нужно составить систему линейных уравнений (3.1.3). Ее решение позволит определить значения коэффициентов полинома, удовлетворяющего условиям (3.1.2). Из единственности решения системы (3.1.3) следует единственность интерполяционного полинома Pn(x). Зная интерполяционный полином, можно определить его значение Pn(x*) в любой, не совпадающей с узлами сетки hx точке x* Î (a, b), осуществив тем самым интерполяцию таблично заданной функции f(x) в точку x*.

Рассмотренный способ построения интерполяционного полинома
с практической точки зрения не всегда эффективен. Существуют альтернативные способы его построения, не требующие решения системы линейных уравнений (3.1.3). Эти способы более подробно рассмотрены в § 3.2 и 3.3.

Существенным недостатком полиномиальной интерполяции является то, что при большом числе узлов сетки степень полинома Pn(x), определяемого на отрезке [a, b], будет большой. Известно, что полиномы с большими степенями крайне неустойчиво ведут себя в промежутках между узлами сетки hx, часто принимая там неоправданно большие по абсолютной величине значения, заведомо не соответствующие возможным значениям функции f(x). Поэтому на практике не рекомендуется строить интерполяционные полиномы выше 5-6-й степени.

В связи с этим при большом количестве узлов сетки hx целесообразно
либо строить интерполяционный полином для небольшой части узлов сетки, находящихся на оси Ох с обеих сторон от точки, в которую предполагается интерполировать данные, либо применять неполиномиальную интерполяцию таблично заданных функций. Широкое распространение получили методы интерполяции так называемыми сплайнами. Пример постановки и решения задачи сплайн-интерполяции приведен в § 3.4.

Рассмотрим другую формальную постановку задачи аппроксимации таблично заданной функции функцией из заданного класса. Пусть, как и ранее, на отрезке [a, b] задана одномерная сетка

hx = {xi / xi = x– 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).Из заданного класса функций требуется найти такую функцию j (x, с) (где c Î Еm – вектор определяемых параметров), чтобы суммарное отклонение значений функции j (x) в узлах сетки hx от заданных в этих узлах значений функции f(x) было минимальным. Для оценки суммарного отклонения может быть использован, например, критерий метода наименьших квадратов:

F(j (x, c)) = [f(xk) – j (xk, c)]2 Þ min.                                      (3.1.4)

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

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

Более подробно задача аппроксимации таблично заданной функции методом наименьших квадратов рассматривается в § 3.5.

§ 3.2. Интерполяционный полином в форме Лагранжа

Пусть на отрезке [a, b] задана одномерная сетка

hx = {xi / xi = x– 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).

Будем искать интерполяционный полином

Pn(x) = a0 + a1x + a2x2 + … + anxn                                                      (3.2.1)

в виде

                                                                       (3.2.2)

Для выполнения условий интерполяции по узлам сетки hx

Pn(xi) = yi, i = 0, 1, 2, …, n                                                                 (3.2.3)

достаточно, чтобы

                         (3.2.4)

Действительно, если функции Ck(x), k = 0, 1, 2, …, n удовлетворяют условиям (3.2.4), то и полином Pn(x) в виде (3.2.2) будет удовлетворять условиям интерполяции (3.2.3).

Так как интерполяционный полином Pn(x) является полиномом степени n, то с учетом (3.2.2) функции Ck(x), k = 0, 1, 2, .., n можно задавать полиномами той же степени.

С учетом условий (3.2.4) будем определять Ck(x), k = 0, 1, 2, …, n
в виде

  (3.2.5)

Очевидно, что при таком задании функций Ck(x), k = 0, 1, 2, …, n
условия (3.2.4) будут выполнены.

Подставив (3.2.5) в (3.2.2), получим следующее выражение для Pn(x), x Î (a, b):

   (3.2.6)

Введем в рассмотрение функцию w(x):

w(x) = (xx0)(xx1)(xx2)(xx3) … (xxn).                                (3.2.7)

Очевидно, что w(x) является полиномом степени n + 1. Продифференцировав w(x), найдем выражение для w¢(x) в узлах сетки hx:

w¢(xk) = (xkx0)(xkx1) … (xkxk – 1)(xkxk + 1) … (xkxn),

k = 0, 1, 2, …, n.                                                                                 (3.2.8)

Подставив (3.2.7) и (3.2.8) в (3.2.6), перепишем выражение для интерполяционного полинома в виде

                                                        3.2.9)

Полином в виде (3.2.9) называют интерполяционным полиномом в форме Лагранжа.

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

                           (3.2.10)

где xk, k = 0, 1, 2, …, n – узлы сетки hx.

Произведение элементов, находящихся на главной диагонали матрицы Аh(x), позволяет вычислить значение w(x) в точке x, а произведение элементов k-й строки (исключая диагональный элемент) этой матрицы – значения w¢(xk), k = 0, 1, 2, …, n. В результате достаточно легко найти значения коэффициентов

                                       (3.2.11)

для интерполяционного полинома в форме Лагранжа.

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

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

hx = {xi  / xi = x– 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

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

Задать x* Î (a, b).

1. Положить k: = 0, S: = 0.

2. Вычислить элементы матрицы Аh(x) для заданной сетки hx и x = x*.

3. Вычислить значение w(x*) как произведение диагональных элементов матрицы Аh(x).

4. Вычислить значение w¢(xk), как произведение всех элементов k-й строки матрицы Аh(x), за исключением диагонального элемента.

Рассчитать значение коэффициента

5. Рассчитать r = Ck(x*) f(xk).

6. Положить S: = S + r.

7. Если k = n, перейти к выполнению п. 9.

8. Положить k: = k + 1, перейти к выполнению п. 4.

9. Процесс завершен. Полученное значение S и есть результат интерполяции табличных данных в точку x*.

Для визуальной оценки качества интерполирования желательно вычислить значения интерполяционного полинома Pn(x) с достаточно малым шагом по переменной x и построить его график на отрезке [a, b], что позволит наглядно сопоставить его с табличными данными.

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

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

2. Как уже было сказано ранее, не рекомендуется использовать интерполяционные полиномы выше 6-го порядка: в этом случае при их построении может быть использовано в совокупности не более семи узлов сетки, находящихся на оси Ох с обеих сторон от точки x* Î (a, b).

3. Имеет место следующая оценка точности полиномиальной интерполяции:

                                                            (3.2.12)

где Mn + 1 – верхняя граница значения модуля (n + 1)-й производной функции f(x) на отрезке [a, b].

В соответствии с (3.2.12) погрешность интерполяции будет равна нулю, если интерполируемая функция представляет собой полином, степень которого не превышает n.

4. Интерполяционный полином в форме Лагранжа часто оказывается удобным для проведения различных теоретических исследований в области вычислительной математики.

§ 3.3. Интерполяционный полином в форме Ньютона

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

hx = {xi  / xi = x– 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы значения yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).

Интерполяционный полином в форме Ньютона для x Î (a, b) может быть записан в виде

Pn(x) = f(x0) + f(x0, x1)×(xx0) + f(x0, x1, x2)×(xx0)(xx1) + …

…+ f(x0, x1, x2, …, xn – 1, xn)×(xx0)(xx1) … (xxn - 2)(xxn - 1)    (3.3.1)

или в более компактной форме

,                                 (3.3.2)

где f(x0, x1,…, xk) – разделенная разность k-го порядка.

Разделенные разности для табличных данных рассчитываются следующим образом:

·          разделенные разности 1-го порядка по формуле

·          разделенные разности 2-го порядка через разделенные разности 1-го порядка:

В общем случае разделенные разности k-го порядка определяются
через разделенные разности (k – 1)-го порядка:

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

Нетрудно показать, что для интерполяционного полинома Pn(x)
в форме Ньютона выполняются условия интерполяции по узлам сетки:

Pn(xk) = f(xk), k = 0, 2, …, n.

Кроме того, можно убедиться в том, что интерполяционный полином в форме Ньютона только по форме отличается от интерполяционного полинома в форме Лагранжа, представляя собой на заданной сетке один и тот же интерполяционный полином. Поскольку, как было отмечено ранее, интерполяционный полином единственен, то последнее утверждение является вполне естественным и должно выполняться для любых форм представления этого полинома, если сетка и значения функции не изменяются.

В соответствии с соотношением (3.3.1) процесс интерполяции таблично заданной на отрезке [a, b] функции в заданную точку x* Î (a, b)
с помощью интерполяционного полинома в форме Ньютона можно представить в виде следующего алгоритма:

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

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

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

Задать x* Î (a, b).

1. Положить S: = f(x0); D: = 1; k: = 1.

2. Вычислить разделенные разности k-го порядка

f(xi, xi + 1, xi + 2, …, xi + k) для всех i = 0, 1, 2, …, nk;

D: = D×(x*xk – 1).

3. Вычислить r = f(x0, x1, x2, …, xk)×D.

4. Вычислить S: = S + r.

5. Положить k: = k + 1.

6. Если k £ n, перейти к выполнению п. 2.

7. Процесс завершен. S = Pn(x*) – значение интерполяционного полинома в точке x*.

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

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

Таблица 3.3.1

№ узла

xi

f(xi)

f(xi, xi+1)

f(xi, xi+1, xi+2)

f(xi, xi+1, xi+2, xi+3)

f(xi, xi+1, xi+2, …, xi+n)

0

x0

f(x0)

 

 

 

 

 

 

 

f(x0, x1)

 

 

 

1

x1

f(x1)

 

f(x0, x1, x2)

 

 

 

 

 

f(x1, x2)

 

f(x0, x1, x2, x3)

 

2

x2

f(x2)

 

f(x1, x2, x3)

 

 

 

 

 

f(x2, x3)

 

f(x1, x2, x3, x4)

 

3

x3

f(x3)

 

f(x2, x3, x4)

 

 

 

 

 

f(x3, x4)

 

f(x2, x3, x4, x5)

 

4

x4

f(x4)

 

f(x3, x4, x5)

 

f(x0, x1, x2,…, x8)

 

 

 

f(x4, x5)

 

f(x3, x4, x5, x6)

 

5

x5

f(x5)

 

f(x4, x5, x6)

 

 

 

 

 

f(x5, x6)

 

f(x4, x5, x6, x7)

 

6

x6

f(x6)

 

f(x5, x6, x7)

 

 

 

 

 

f(x6, x7)

 

f(x5, x6, x7, x8)

 

7

x7

f(x7)

 

f(x6, x7, x8)

 

 

 

 

 

f(x7, x8)

 

 

 

8

x8

f(x8)

 

 

 

 

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

2. Использование интерполяционного полинома в форме Ньютона оказывается удобным, например, когда появляются дополнительные узлы. Так, если добавлен один новый узел, то в формуле (3.3.1) для интерполяционного полинома в форме Ньютона добавится еще одно слагаемое, которое нужно прибавить к уже рассчитанной сумме предыдущих слагаемых (они, очевидно, не изменяются). Для расчета значения дополнительного слагаемого в нижней части таблицы для конечных разностей добавляется еще одна строка, соответствующая этому узлу. Ранее рассчитанные разности пересчитывать не нужно, необходимо лишь определить дополнительно значения тех разделенных разностей, которые будут находиться на нижней стороне нового "треугольника" разделенных разностей.

3. Формулу (3.3.1) для интерполяционного полинома в форме Ньютона можно записать для любого способа упорядочивания узлов. Так, например, при упорядочивании узлов в обратном порядке (xn, xn – 1, xn – 2, …, x2, x0) формула (3.3.1) для представления интерполяционного полинома в форме Ньютона примет вид

Pn(x) = f(xn) + f(xn, xn – 1)×(xxn) + f(xn, xn – 1, xn – 2)×(xxn)(xxn – 1)+…

…+ f(xn, xn – 1, xn – 2,… x1, x0)×(xxn)(xxn – 1)… (xx2)(xx1).     (3.3.3)

Если интерполяционный полином в форме Ньютона Pn(x) определяется по формуле (3.3.1), его называют интерполяционным полиномом для интерполирования вперед (или для интерполирования в начале таблицы). Этот полином целесообразно использовать, когда точка x* находится в начале таблицы.

Если Pn(x) определяется по формуле (3.3.3), его называют интерполяционным полиномом для интерполирования назад (или для интерполирования в конце таблицы). Его целесообразно использовать, когда точка x* находится в нижней части таблицы.

4. Как уже говорилось, не рекомендуется использовать интерполяционные полиномы выше 6-го порядка: в этом случае при их построении может быть использовано в совокупности не более семи узлов сетки, находящихся на оси Ох с обеих сторон от точки x* Î (a, b).

5. Оценка точности интерполяции с помощью интерполяционного полинома в форме Ньютона совпадает с оценкой точности интерполяции, приведенной для интерполяционного полинома в форме Лагранжа, поскольку оба этих полинома представляют собой один и тот же интерполяционный полином и различаются только формой записи:

                                                         (3.3.4)

где Mn + 1 – верхняя граница значения модуля (n + 1)-й производной функции f(x) на отрезке [a, b].

В соответствии с (3.3.4) погрешность интерполяции будет равна нулю, если интерполируемая функция представляет собой полином, степень которого не превышает n.

§ 3.4. Сплайн-интерполяция таблично заданной функции

Интерполяция таблично заданных функций с помощью сплайнов "лишена" главного недостатка полиномиальной интерполяции – роста степени интерполяционного полинома при увеличении количества узлов сетки. Сплайн-интерполяция представляет собой кусочно-полиномиальную интерполяцию, которая во многих случаях оказывается более эффективной, чем попытки подобрать один интерполяционный полином для отрезка [a, b], на котором представлены табличные данные. Следует отметить, что идеи использования кусочно-полиномиального приближения функций в вычислительной математике возникли достаточно давно. Классическим примером можно считать идеи, положенные в основу метода Эйлера для решения задачи Коши для обыкновенных дифференциальных уравнений первого порядка. Современная теория сплайн-интерполяции и прикладной математический аппарат для использования сплайнов с целью интерполяции и аппроксимации функций начали интенсивно развиваться с середины ХХ веке. В настоящее время это общепризнанное и все еще развивающееся направление вычислительной математики с достаточно широкой областью практического применения. При решении задач интерполяции преимущество сплайнов перед интерполяционными полиномами заключается в гарантированных устойчивости вычислительного процесса и сходимости сплайн-интерполяции при увеличении количества узлов сетки. В настоящем пособии рассмотрена задача интерполирования таблично заданной на отрезке [a, b] функции f(x) с использованием кубических сплайнов, которые широко применяются для интерполяции функций.

Пусть на отрезке [a, b] задана одномерная сетка

hx = {xi  / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы значения yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x). Очевидно, что все узлы сетки различны и упорядочены по возрастанию.

Кубическим сплайном будем называть функцию S(x), заданную
на отрезке [a, b] и удовлетворяющую следующим условиям:

1. S(xi) = yi, i = 0, 1, 2, …, n.

2. S(x), S ¢(x) и S ¢¢(x) непрерывны на отрезке [a, b].

На любом отрезке [xi – 1, xi], i = 1, 2, …, n  S(x) представляет собой полином 3-й степени:

Si(x) = ai + bi(xxi) + (ci / 2)(xxi)2 + (di / 6)(xxi)3,

x Î [xi – 1, xi], i = 1, 2, …, n.                                                                (3.4.1)

Для построения кубического сплайна на отрезке [a, b] необходимо найти коэффициенты ai, bi, ci, di, i = 1, 2, …, n – всего 4n неизвестные величины. Для их нахождения в соответствии с определением кубического сплайна имеются условия:

– совпадения значений S(x) и таблично заданной функции f(x) в узлах сетки hx:

S(xi) = yi , i = 0, 1, 2, …, n;

– непрерывности функции S(x), ее первой и второй производных.

Эти условия сводятся к необходимости обеспечения непрерывности S(x) и ее производных во всех внутренних узлах (точках) xi, i = 1, 2, …,
n – 1 сетки hx, поскольку в этих точках происходит смена аналитического задания кусочно-полиномиальной функции S(x), в остальных точках отрезка [a, b] указанные функции непрерывны по определению.

Таким образом, для нахождения 4n неизвестных коэффициентов ai, bi, ci, di, i = 1, 2, …, n имеется только (n + 1) + 3(n – 1) = 4n – 2 условий. Два недостающих условия получаются из так называемых граничных условий, которые назначаются исходя из физических или других соображений, связанных с особенностями интерпретации таблично заданной функции f(x). Выберем в качестве граничных условий равенство нулю второй производной функции S(x) в граничных точках отрезка [a, b]:

S¢¢(a) = S¢¢(b) = 0.

В результате получим систему линейных уравнений из 4n уравнений для определения 4n неизвестных коэффициентов ai, bi, ci, di, i = 1, 2, …, n. Предварительный анализ этих уравнений и ряд несложных преобразований приводят к достаточно простой последовательности операций для нахождения значений указанных коэффициентов.

Коэффициенты ai, i = 1, 2, …, n определяются из соотношений

ai = yj, i = 1, 2, …, n.                                                                          (3.4.2)

Для нахождения коэффициентов ci, i = 1, 2, …, n необходимо решить трехдиагональную систему линейных уравнений

 3.4.3)

Очевидно, что эта система имеет матрицу с диагональным преобладанием, у нее единственное решение, для нахождения которого может быть использован метод прогонки.

После того как коэффициенты ci, i = 1, 2, …, n будут найдены, коэффициенты di, i = 1, 2, …, n могут быть определены по формуле

di = (ci ci  1) / hi, i = 1, 2, …                                                             (3.4.4)

Наконец, находим коэффициенты bi, i = 1, 2, …, n по формуле

                               (3.4.5)

В результате вычислений будет построен интерполяционный кубический сплайн S(x).

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

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

hx = {xi / xi = xi –1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

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

Задать x* Î (a, b).

1. Положить ai = yj, i = 0, 1, 2, …, n.

3. Составить и решить трехдиагональную систему (3.4.3) методом прогонки. Определить значения коэффициентов ci, i = 0, 1, 2, …, n.

4. Определить значения коэффициентов di и bi, i = 1, 2, 3, …, n, воспользовавшись формулами (3.4.4) и (3.4.5) соответственно.

5. Определить значение индекса 0 < k £ n из условия x* Î [xk – 1, xk].

6. Вычислить по формуле

S(x*) = Sk(x*) = ak + bk(x*xk) + (ck / 2)(x*xk)2 + (dk / 6)(x*xk)3.

7. Процесс завершен: S(x*) – результат интерполяции табличных данных в точку x* Î (a, b).

Результаты расчета коэффициентов кубического сплайна удобно представлять в виде таблицы, аналогичной табл. 3.4.1.

 

 

 

 

Таблица 3.4.1

k

xk

yk

ak

bk

Ck сk

dk

0

x0

y0

a0

-

0

-

1

x1

y1

a1

b1

C1

d1

2

x2

y2

a2

b2

C2

d2

3

x3

y3

a3

b3

C3

d3

n – 1

xn – 1

yn – 1

an – 1

bn – 1

cn – 1

dn – 1

n

xn

yn

an

bn

0

dn

Для визуальной оценки качества интерполирования желательно вычислить значения сплайна S(x) с достаточно малым шагом по переменной x и построить его график на отрезке [a, b], сопоставив с аналогичным графиком интерполяционного полинома Pn(x) и данными таблицы.

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

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

                                                         (3.4.6)

где M4 – максимум модуля четвертой производной интерполируемой функции на отрезке [a, b].

Предполагается, что на отрезке [a, b] функция f(x) непрерывна вместе со своими производными до 4-го порядка включительно. В соответствии с (3.4.6) погрешность интерполяции будет равна нулю, если интерполируемая функция представляет собой полином, степень которого
не превышает 3.

2. Интерполирование кубическими сплайнами является сходящейся процедурой в том смысле, что при стремлении максимального шага сетки к нулю (при этом, очевидно, неограниченно возрастает количество узлов сетки) погрешность сплайн-интерполяции равномерно стремится к нулю на отрезке [a, b] не только для интерполируемой функции f(x), но и для ее первых двух производных.

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

4. С увеличением количества узлов сетки hx на отрезке [a, b] степень используемых при сплайн-интерполяции полиномов не изменяется. Прямо пропорционально количеству узлов увеличивается только количество составляющих сплайн S(x) кубических полиномов Sk(x), k = 1, 2, 3, …, n. Как следствие, растет порядок решаемой для определения коэффициентов ci,
i = 0, 1, 2, …, n трехдиагональной системы линейных уравнений и общее количество вычислительных операций, требуемых для интерполирования табличных данных в заданную точку.

5. Как видно из описания построения кубического сплайна, он
не может быть однозначно определен из условий интерполяции; поэтому необходимо ввести два дополнительных условия, что создает определенные сложности при практической интерполяции. В любом случае при использовании стандартных компьютерных процедур сплайн-интерполяции желательно представлять, какие именно дополнительные условия в них используются. Кубический сплайн, при построении которого дополнительно вводятся условия

S ¢¢(a) = S ¢¢(b) = 0,

часто называют естественным.

6. Интерполирование с помощью сплайнов широко используется на практике. Однако нельзя сказать, что фактические результаты интерполяции всегда бесспорно удовлетворяют исследователей. Как отмечается
в литературе, имеется много случаев, когда сплайн не вполне удовлетворительно позволяет восстанавливать функции, для которых, например,
характерны быстрые изменения значений на малых промежутках. С целью устранения подобных недостатков разрабатываются различные модификации классических процедур, с которыми можно ознакомиться в соответствующей литературе.