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

Интерполяция таблично заданных функций с помощью сплайнов "лишена" главного недостатка полиномиальной интерполяции – роста степени интерполяционного полинома при увеличении количества узлов сетки. Сплайн-интерполяция представляет собой кусочно-полиномиальную интерполяцию, которая во многих случаях оказывается более эффективной, чем попытки подобрать один интерполяционный полином для отрезка [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
x0 y0 a0 - -
x1 y1 a1 b1 C1 d1
x2 y2 a2 b2 C2 d2
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 dn

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

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

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

(3.4.6)

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

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

return false">ссылка скрыта

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

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

Представим ситуацию, когда табличные данные имеют заметную вычислительную погрешность. В других случаях возникает необходимость подбора на основе имеющихся экспериментальных данных, представленных в табличной форме, некоторой функциональной зависимости из заданного класса функций, свойства которого наиболее соответствуют предполагаемым физическим свойствам изучаемых процессов или позволяют выявить характерную для этих процессов тенденцию. С математической точки зрения в этих и других возможных аналогичных случаях имеется функция с неизвестным (или известным, но слишком сложным с точки зрения проводимого исследования) аналитическим заданием, представленная в узлах заданной сетки значениями, которые были получены в результате прямых или косвенных измерений или некоторого вычислительного эксперимента. Такие данные будем называть экспериментальными. Поскольку по ряду причин (например, в силу большой погрешности) интерполирование такой таблично заданной функции не является целесообразным, ставится задача ее аппроксимации (замены) функцией из заданного класса, наиболее согласованной в смысле некоторого критерия с имеющимися экспериментальными данными. Подбираемую функцию будем называть моделью, а ее значения в узлах заданной сетки – модельными данными. Если в качестве критерия согласованности выбирается критерий минимума суммы квадратов отклонений модельных данных от экспериментальных, говорят об аппроксимации таблично заданной функции методом наименьших квадратов (МНК).

Формальная постановка задачи осуществляется следующим образом. Пусть на отрезке [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). Пусть также для аппроксимации табличных данных выбран некоторый класс функций j (х, c0, c1, c2, …, cm), m < n, где c0, c1, c2, …, cm – неопределенные коэффициенты, выбор значений которых позволяет определить конкретную функцию из выбранного класса. Требуется найти значения коэффициентов c0, c1, c2,…, cm, для которых выполнено условие j (хi, c0, c1, c2, …, cm))2 Þ min. (3.5.1)

Критерий (3.5.1), на основании которого осуществляется выбор значений коэффициентов c0, c1, c2, …, cm, является стандартным, используемым в МНК. Выбранные в соответствии с этим критерием значения коэффициентов позволяют определить среди множества функций конкретную функцию, наиболее согласованную с табличными (экспериментальными) данными или, иначе, обеспечивающую наилучшее среднеквадратическое приближение.

Как отмечалось выше, функция j (х, c0, c1, c2, …, cm) называется моделью, тогда коэффициенты c0, c1, c2,…, cm будем называть параметрами
модели или модельными. В дальнейшем ограничимся рассмотрением случая, когда модель линейно зависит от параметров и ее можно представить в виде

j (х, c0, c1, c2, …, cm) = c0j0(х) + c1j1(х) + c2j2(х) + … + cmjm(х). (3.5.2)

Модель вида (3.5.2) часто называют обобщенным полиномом. Здесь j0(х), j1(х), j2(х), jm(х) – множество так называемых базисных функций. Базисные функции могут быть как линейными, так и нелинейными функциями переменной x. Независимо от этого модель (3.5.2) остается линейной, поскольку она линейно зависит от модельных параметров c0, c1, c2, …, cm.
В качестве базисных могут быть выбраны, например, степенные функции

j0(х) = 1, j1(х) = х, j2(х) = х2, …, jm(х) = хm.

Тогда модель будет представлять собой полином степени m:

j (х, c0, c1, c2, …, cm) = c0 + c1х + c2х2 + … + cmхm. (3.5.3)

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

Итак, для линейной модели (3.5.2) требуется найти значения модельных параметров c0, c1, c2, …, cm, обеспечивающих выполнение условия (3.5.1). Эту задачу можно решить двумя способами.