40. Интерполяционные сплайны. Общие понятия. Постановка задачи. Кубические сплайны.

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

Интерполяция таблично заданных функций с помощью сплайнов "лишена" главного недостатка полиномиальной интерполяции – роста степени интерполяционного полинома при увеличении количества узлов сетки. Сплайн-интерполяция представляет собой кусочно-полиномиальную интерполяцию, которая во многих случаях оказывается более эффективной, чем попытки подобрать один интерполяционный полином для отрезка [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. Интерполирование с помощью сплайнов широко используется на практике. Однако нельзя сказать, что фактические результаты интерполяции всегда бесспорно удовлетворяют исследователей. Как отмечается
в литературе, имеется много случаев, когда сплайн не вполне удовлетворительно позволяет восстанавливать функции, для которых, например, характерны быстрые изменения значений на малых промежутках. С целью устранения подобных недостатков разрабатываются различные модификации классических процедур, с которыми можно ознакомиться в соответствующей литературе.