MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Spline Toolbox

Построение интерполяционных сплайнов в B-форме

Интерполяция сплайном в B-форме позволяет просто управлять гладкостью сплайна. Сплайн в B-форме определяется последовательностью узлов {ti}, среди которых могут быть повторяющиеся. Число повторений узла определяет гладкость сплайна в точках разрыва . Во внутренних точках разрыва должно выполняться следующее правило:

число повторений узлов в последовательности + число условий непрерывности в точке разрыва = порядку сплайна.

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

Интерполирование при помощи функции spapi

Рассмотрим простой пример. Проинтерполируем функцию x sin x на отрезке [0, 3] сплайном f(x) четвертого порядка (k=4), выбрав семь полюсов интерполяции

.

Обычно первые k узлов совпадают с первым полюсом, а последние k узлов - с последним. Возьмем следующую последовательность из одиннадцати узлов (число узлов должно равняться числу полюсов интерполяции + порядок сплайна, в данном случае 11=7+4)

.

и применим функцию spapi:
tau = [0  0.5  1.1   1.8   2.2    2.8    3]; % полюса интерполяции
y = tau.*sin(tau); % значения интерполируемой функции в полюсах интерполяции
t = [0 0 0 0    1.2     1.7    2.3    3 3 3 3];  % последовательность узлов сплайна 
sp = spapi(t, tau, y); % конструирование сплайна
figure 
fnplt(sp)  % построение графика сплайна
hold on
plot(tau, y, 'ro') % отображение данных маркерами
Получаем следующий график (см. рис. 1), где сплайн f(x) нарисован линией, а маркеры соответствуют значениям интерполируемой функции x sin x в полюсах интерполяции :


Рис. 1. Интерполяция сплайном f(x) в B-форме при помощи spapi.

Убедимся, что внутренние узлы являются точками разрыва. Их кратность равна единице, следовательно, число условий непрерывности (порядок сплайна минус кратность узла, т.е 4-1) в каждой из них равно трем и, следовательно, должна быть непрерывна сама сплайн-функция, ее первая и вторая производные. Для этого построим графики первой, второй и третьей производных сплайна f(x), вычислив значения производных сплайна при помощи функции fnder. На графики поместим вертикальные линии в точках разрыва (см. рис. 2):
figure
subplot(3,1,1)
fnplt(fnder(sp,1)) % находим первую производную сплайна и строим ее график
title('1st derivative of the spline')
hold on
lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей
plot([t(5:7); t(5:7)], [ymin ymax],'k') % рисуем линии в точках разрыва
subplot(3,1,2)
fnplt(fnder(sp,2)) % находим вторую производную сплайна и строим ее график
title('2nd derivative of the spline')
hold on
lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей
plot([t(5:7); t(5:7)], [ymin ymax],'k') % рисуем линии в точках разрыва
subplot(3,1,3)
fnplt(fnder(sp,3)) % находим третью производную сплайна и строим ее график
title('3rd derivative of the spline') % рисуем линии в точках разрыва
hold on
lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей
plot([t(5:7); t(5:7)], [ymin ymax],'k') 

Рис. 2. Графики первой, второй и третьей производной сплайна f(x) в B-форме.

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

.

t = [0 0 0 0   1.2    1.7    1.7     3 3 3 3]
sp = spapi(t, tau, y);
Графики сплайна и его производных строится аналогично, из графиков производных видно (см. рис. 3), что вторая производная сплайна имеет разрыв в узле 1.7:


Рис. 3. Графики первой, второй и третьей производных сплайна в B-форме для другой последовательности узлов (узел 1.7 имеет кратность 2).

B-форма, базисные сплайны

Сплайн в B-форме является суммой B-сплайнов (базисных сплайнов), каждый из которых отличен от нуля на некотором небольшом интервале. Строгое определение B-сплайнов и их свойства мы рассмотрим позже. Сейчас важно понять идею, заложенную в основу B-формы и B-сплайна. Вернемся к нашему примеру интерполяции функции на отрезке [0, 3] кубическим сплайном (т.е. сплайном 4-го порядка) в B-форме с непрерывными 1-ой и 2-ой производными в каждой внутренней точке с узлами


Сплайн f(x) должен совпадать с функцией функции x sin x в полюсах интерполяции


tau = [0  0.5  1.1   1.8   2.2    2.8    3]; % полюса интерполяции
y = tau.*sin(tau); % значения интерполируемой функции в полюсах интерполяции
t = [0 0 0 0    1.2     1.7    2.3    3 3 3 3];  % последовательность узлов сплайна 
sp = spapi(t, tau, y) % конструирование сплайна
Обратимся к содержимому структуры sp, в которую записана информация о B-форме сплайна f(x):
sp = 
      form: 'B-'
     knots: [0 0 0 0 1.2000 1.7000 2.3000 3 3 3 3]
     coefs: [4.2538e-019 -0.0368 0.8292 1.8139 1.9109 1.0881 0.4234]
    number: 7
     order: 4
       dim: 1
Самыми важными для нас сейчас является поля number и coefs - оказывается B-форма состоит из семи сплайнов (sp.number=7), назовем их B1(x), B2(x), ... , B7(x). Сплайн-функция f(x) является суммой этих сплайнов, умноженных на соответствующие коэффициенты Ci:


Массив коэффициентов Ci записан в sp.coefs. Что же это за сплайны B1(x), B2(x), ... , B7(x)? Оказывается каждый из них построен по своей последовательности из пяти узлов, указанной в таблице:


т.е. сплайн Bi(x) построен по узлам ti, ti+1, ..., ti+k. Определение B-сплайна через его узлы мы рассмотрим позже. Сейчас нам понадобится функция spmak, входящая в состав Spline Toolbox, которая конструирует B-сплайн по заданной последовательности узлов. Для построения, например, сплайна Bi(x) следует обратиться к ней так:
B1 = spmak([t(1) t(2) t(3) t(4) t(5)], 1),
или задать узлы при помощи сечения массива t:
  
B1 = spmak([t(1:5)], 1).
Результатом является структура B1 с информацией о сплайне. Создадим сплайны B1(x), B2(x), ... , B7(x), запишем их в структуры B1, B2, ..., B7, соответственно, и построим их графики разыми цветами (см. рис. 4) при помощи функции fnplt (отдельные структуры для сплайнов взяты только для наглядности изложения, можно было ввести массив структур и все сделать в циклах):
% Создание сплайнов B1, B2, ... , B7
B1 = spmak(t(1:5), 1);    B2 = spmak(t(2:6), 1);    B3 = spmak(t(3:7), 1);
B4 = spmak(t(4:8), 1);    B5 = spmak(t(5:9), 1);    B6 = spmak(t(6:10), 1);
B7 = spmak(t(7:11), 1);
% Построение графиков сплайнов B1, B2, ... , B7
figure
fnplt(B1, 'r'); hold on; 
fnplt(B2, 'g'); fnplt(B3, 'b'); fnplt(B4, 'c');
fnplt(B5, 'm'); fnplt(B6, 'y');  fnplt(B7, 'k');
% Подписываем координаты узлов сплайна
set(gca, 'XTick', t(4:8), 'XTickLabel', num2str(t(4:8)'), 'XGrid', 'on')

% Добавление легенды
str = {'spline B_1(x)'; 'spline B_2(x)'; 'spline B_3(x)';...
    'spline B_4(x)'; 'spline B_5(x)'; 'spline B_6(x)'; 'spline B_7(x)'};
legend(str, -1)

Рис. 4. Графики B-сплайнов, образующих сплайн f(x).

Вне интервала (ti, ti+k) B-сплайн Bi(x) тождественно равен нулю, в чем несложно убедиться, построив их графики (см. рис. 5), например:
figure
% строим график B1 на отрезке [-1, 4]
subplot(2, 1 ,1)
fnplt(B1, [-1 4], 'r') 
title('spline B_1')
set(gca,'XTick', [t(1) t(5)], 'XTickLabel', num2str([t(1) t(5)]'), 'XGrid', 'on')
% строим график B2 на отрезке [-1, 4]
subplot(2, 1, 2)
fnplt(B2, [-1 4], 'g') 
title('spline B_2')
set(gca, 'XTick', [t(2) t(6)], 'XTickLabel', num2str([t(2) t(6)]'), 'XGrid', 'on')

Рис. 5. Графики сплайнов B1(x) и B2(x) на отрезке [-1, 4].

Обратите внимание, что сплайн B1(x) имеет разрыв в точке 0 потому, что его среди его узлов t1=0, t2=0, t3=0, t4=0, t5=1.2 узел со значением 0 встречается 4 раза. Так как сам сплайн 4-го порядка, то число условий непрерывности равно нулю (4-4=0), т.е. нет непрерывности не только производных, но и самого сплайна.
Сумма сплайнов B1(x), B2(x), ... , B7(x), умноженных на соответствующие коэффициенты (которые возвращаются в sp.coefs при вызове sp = spapi(t, tau, y)),


является сплайном в B-форме, ровно тем, который возвращает функция spapi. Это несложно проверить, вычислив данную сумму сплайнов. Для нахождения значений каждого из сплайнов B1(x), B2(x), ... , B7(x) в точках отрезка [0, 3] с шагом 0.05 применим функцию fnval, значения суммы в этих точках запишем в вектор f и построим график зависимости суммы от x на отрезке [0, 3]:
% записываем коэффициенты сплайна в B-форме в вектор c
c = sp.coefs;
% вычисляем значения сплайна на отрезке [0 3] с шагом 0.05 
x = 0:0.05:3;
f = c(1)*fnval(B1, x) + c(2)*fnval(B2, x) + c(3)*fnval(B3, x) + c(4)*fnval(B4, x) + ...
    c(5)*fnval(B5, x)+c(6)*fnval(B6, x) + c(7)*fnval(B7, x);
% строим графики сплайна и исходных данных
figure
plot(x, f)
hold on
plot(tau, y, 'ro')
Построенный график совпадает с графиком сплайна, полученного ранее при непосредственном использовании функции spapi для интерполяции функции x sin x, поэтому график сплайна здесь не приводится (см. рис. 1).

Выбор узлов B-сплайнов

Как мы уже замечали, последовательность узлов B-формы сплайна не может выбираться независимо от последовательности полюсов интерполяции. Во-первых, число узлов должно равняться числу полюсов + порядок сплайна. Во-вторых, для существования и единственности сплайна требуется выполнение еще условия теоремы Шенберга-Уитни.
Начнем с примера. Попробуем взять такую последовательность узлов:


а полюса интерполяции оставим теми же самыми


При выполнении команд
tau = [ 0   0.5  1.1  1.8  2.2  2.8      3];
y = tau.*sin(tau);
t = [0 0 0 0     0.3     0.4    1.9    3 3 3 3]
sp = spapi(t, tau, y);
получим ошибку и сообщение о том, что матрица является вырожденной. Дело в том, что при интерполяции B-сплайнами для нахождения коэффициентов B-формы необходимо решить некоторую систему линейных алгебраических уравнений, матрица которой в данном случае оказывается вырожденной. Для однозначной разрешимости этой системы линейных уравнений, а следовательно, и для однозначного определения интерполяционного сплайна все элементы последовательности полюсов интерполяции и узлов должны удовлетворять неравенствам:


в которых строгое равенство допускается лишь для крайних узлов (разумеется, имеется ввиду, что последовательность узлов упорядочена по возрастанию).
В данном случае уже не выполняется условие , так как

Конечно, для построения подходящей последовательности узлов сплайна не требуется подбирать их вручную. В Spline Toolbox входит функция aptknt, которая по заданной последовательности полюсов интерполяции tau возвращает последовательность t узлов, обеспечивающую существование и единственность сплайна заданного порядка k. Обращение к функции aptknt выглядит следующим образом:
t = aptknt(tau, k)
Для интерполяции B-сплайнами 4-го порядка функции в нашем примере подошла бы такая последовательность команд:
tau = [ 0   0.5 1.1 1.8 2.2 2.8      3];
y = tau.*sin(tau);
t = aptknt(tau, 4)
sp = spapi(t, tau, y);
figure 
fnplt(sp)  % построение графика сплайна
hold on
plot(tau, y, 'ro') % отображение данных маркерами
Более того, интерфейс функции spapi позволяет вообще не заботится об узлах сплайна и указать только степень, полюса интерполяции и значения функции в них:
sp = spapi(4,tau,y);
Результат будет тот же самый, поскольку при таком способе обращения к функции spapi она вызывает функцию aptknt для конструирования допустимой последовательности узлов сплайна, удовлетворяющей условиям теоремы Шенберга-Уитни.

Функция aptknt возвращает такую последовательность узлов:


где внутренние узлы сплайна получаются усреднением соответствующих полюсов интерполяции:


для i=1,2,...,n-k. Предполагается, что последовательность полюсов интерполяции содержит не менее элементов, является неубывающей и для всех i. Эти узлы удовлетворяют условиям теоремы Шенберга-Уитни и обеспечивают существование и единственность интерполяционного сплайна.
Вообще говоря, выбор хороших полюсов и узлов интерполяции сплайна (оптимальных в некотором смысле) - отдельная задача. В Spline Toolbox имеется специальная функция optknt, которая так же, как и aptknt, по заданным полюсам интерполяции и заданному порядку сплайна возвращает последовательность его узлов, удовлетворяющую условию теоремы Шенберга-Уитни. Узлы сплайна, которые строит функция optknt, приводят к наилучшему приближению сплайном.
Обратная задача - выбор в некотором смысле оптимальных полюсов для интерполяции сплайном заданного порядка с заданными узлами - решается при помощи функции chbpnt. Хорошее приближение также дают полюса интерполяции, получаемые путем усреднения узлов сплайна. Для сплайна порядка k с последовательностью узлов такие полюса определяются по формуле:


Для их вычисления в Spline Toolbox предназначена функция aveknt. Обсуждению проблемы выбора полюсов интерполяции и узлов сплайна будет посвящен специальный раздел.

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

Функция augknt позволяет просто создать последовательность узлов сплайна по заданным точкам разрыва так, что в каждой точке разрыва узел будет иметь заданную кратность, например
>> t = augknt([-2 -1 0 1 2], 4)
возвращает
t =
    -2    -2    -2    -2    -1     0     1     2     2     2     2
Если требуется задать кратность узла в каждой внутренней точке отдельно, то следует указать вектор кратностей в качестве третьего входного аргумента:
>> t = augknt([-2 -1 1 2], 4, [2 3])

t =
    -2    -2    -2    -2    -1    -1     1     1     1     2     2     2     2
Для конструирования последовательности узлов с заданными кратностями подходит также функция brk2knt, которая требует задания двух входных аргументов - массива точек разрыва и массива их кратностей такой же длины:
>> knots = brk2knt([-1 0 1], [3 2 4])
knots =
    -1    -1    -1     0     0     1     1     1     1
Обратная задача - получение точек разрыва и соответствующих кратностей узлов решается при помощи функции knt2brk:
>>knots = [-1    -1    -1     0     0     1     1     1     1]
>> [breaks,mults] = knt2brk(knots)
breaks =
    -1     0     1
mults =
     3     2     4

Интерполирование B-сплайнами с заданными значениями производных

Функция spapi позволяет конструировать такие B-сплайны, которые принимают заданные значения вместе со своими производными в полюсах интерполяции. Для этого при вызове функции spapi:
sp = spapi(t, tau, y), где t - массив узлов, или
sp = spapi(k, tau, y), где k - порядок сплайна
следует повторить координату полюса в массиве tau и, соответственно, в векторе y задать значения не только сплайна, но и его производной в этом полюсе интерполяции. Для задания значения сплайна и его p-ой производной в некотором полюсе , он должен встречаться p+1 раз в массиве tau, тогда соответствующие значения массива y будут восприниматься как значения сплайна, его первой, второй и т. д., ..., p-ой производной.
Например, для интерполяции кубическим сплайном некоторой функции, заданной таблицей:


подойдет такая последовательность команд (см. результат на рис. 6, на котором приведен график сплайна, его первой производной и маркерами отмечены заданные значения в полюсах интерполяции):
x = [0 1 2 3];
y = [0 -1 0 1];
dy = [2 3 4 2]
sp = spapi(4, [x x], [y dy])
figure
subplot(2, 1, 1)
fnplt(sp)
title('cubic spline')
hold on
plot(x, y, 'ro')
subplot(2,1, 2)
fnplt(fnder(sp))
title('its first derivative')
hold on
plot(x, dy, 'ro')

Рис. 6. Интерполяция B-сплайном по заданным значениям функции и ее производной.

Значения производных интерполируемой функции могут быть заданы не во всех полюсах интерполяции. Рассмотрим пример, в котором требуется интерполировать функцию y=e-xsin3x на отрезке [0, 4] в шести полюсах: так, чтобы значения первой и второй производных сплайна f(x) совпадали со значениями y и в трех полюсах . Для этого задаем массив полюсов интерполяции tau, inline-функцию для вычисления и записываем значения ее первой и второй производных в полюсах в массивы y1 и y2, соответственно:
fun = inline('exp(-x).*sin(3*x)');
tau = 0:5;
y = fun(tau);
y1 = exp(-tau(1:2:5)).*(3*cos(3*tau(1:2:5))-sin(3*tau(1:2:5)));
y2 = -exp(-tau(1:2:5)).*(8*sin(3*tau(1:2:5))+6*cos(3*tau(1:2:5))); 
Далее обращаемся к функции spapi:
f = spapi(4, [tau tau(1:2:5) tau(1:2:5)], [y y1 y2]);
Сравним исходную функцию, полученный сплайн и еще два сплайна: g(x) (который был построен с учетом значений производных ) и сплайн h(x) (который просто совпадает с функцией f(x) в полюсах интерполяции). Результат представлен на рис. 7.
g = spapi(4, [tau tau(1:2:5)], [y y1]);
h = spapi(4, tau, y);

figure
fplot(fun, [0 5], 'b')
hold on
fnplt(f, 'm')
fnplt(g, 'g')
fnplt(h, 'k')
legend('y(x)', 'f(x)', 'g(x)', 'h(x)')

Рис. 7. Интерполяция функции y=e-xsin3x с учетом и без учета значений производных.


Поиск по сайту:

Система Orphus

Яндекс.Метрика