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

Spline Toolbox

Интерполяция кубическими сплайнами, функция csape

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

Здесь находится файл demcsape.m, который содержит все команды этого раздела. Достаточно скопировать его в текущий каталог MATLAB, запустить из командной строки и нажимать <Enter> для получения следующего результата.

Самый простой вариант вызова функции csape выглядит следующим образом:

pp = csape(x, y)

где x - массив со значениями абсцисс, y - массив с соответствующими значениями интерполируемой функции. Функция csape возвращает структуру pp, содержащую информацию о построенном сплайне. Данная структура pp может указываться в качестве входного аргумента других функций Spline Toolbox для визуализации сплайна и выполнения над ним различных операций.

Начнем с интерполяции функции cosx на отрезке в равноотстоящих точках с шагом . Запишем значения независимой переменной в массив x, а функции - в массив y:

x = 0:pi/2:4*pi;
y = cos(x);

Применим теперь функцию csape для интерполяции кубическим сплайном:

pp = csape(x,y)

В командное окно выводится содержимое структуры pp:

pp = 
      form: 'pp'
    breaks: [0 1.5708 3.1416 4.7124 6.2832 7.8540 9.4248 10.9956 12.5664]
     coefs: [8x4 double]
    pieces: 8
     order: 4
       dim: 1

Смысл полей структуры pp мы обсудим чуть позже, а сейчас сравним функцию cosx и построенный сплайн , отобразив их графически на одних осях. Для построения графика сплайна используем функцию fnplt, входящую в Spline Toolbox, а для вывода графика косинуса - функцию fplot:

x = 0:pi/2:4*pi;
y = cos(x);
pp = csape(x, y)
figure
fnplt(pp)
hold on
fplot(@cos, [0 4*pi], 'r')

На график добавим узлы интерполяции и разместим легенду (см. рис. 1):

plot(x, y, 'ok')
legend('cubic spline (pp)', 'cos \itx', 'sites')

Рис. 1. Интерполяция cosx кубическим сплайном (функция csape)

Обсудим теперь смысл кусочно-полиномиальной формы сплайна и ее представление структурой pp, которую возвращает функция csape. Структура pp имеет следующие поля:

·  form - форма спалйна (значение 'pp' говорит о том, что сплайн в кусочно-полиномиальной форме);

·  breaks - вектор точек разрыва (в данном случае 9 точек разрыва [0 1.5708 3.1416 4.7124 6.2832 7.8540 9.4248 10.9956 12.5664], следовательно сплайн составлен из 8 полиномиальных частей);

·  coefs - массив коэффициентов сплайна, каждая строка которого содержит коэффициенты полиномиальных частей (в данном случае массив 8x4, т.к. на каждом из 8-и отрезков полином третьей степени определяется 4-мя коэффициентами);

·  pieces - число полиномов, составляющих сплайн (в данном случае 8);

·  order - порядок сплайна (в данном случае 4, т.к. порядок равен степени плюс 1, т.е. числу коэффициентов полиномов).

·  dim - размерность (в данном случае 1).

Функция MATLAB ployval, предназначенная для вычисления значения полинома по заданному вектору его коэффициентов, позволяет легко построить полиномы, составляющие сплайн. Для сплайна запишем матрицу C с коэффициентами полиномов в массив C, а точки разрыва - в массив ksi:

C = pp.coefs;
ksi = pp.breaks;

На каждом участке для сплайн представляется полиномом третьей степени:

Для наглядности построим графики каждого из полиномов на более широких отрезках , где . Воспользуемся polyval, которая по вектору коэффициентов полинома находит его значения в заданных точках (к примеру, Y=polyval([-3 1 2 3], -1:0.1:1) вычисляет полином в точках от -1 до 1 с шагом 0.1). Кроме того, отметим точки разрыва (см. рис. 2, где каждому участку сопоставлен цвет линии графика соответствующего полинома):

figure
plot(x, y, 'ok')  % выводим точки разрыва маркерами
hold all
plot([x; x], [-1 1], 'k:')  % выводим вертикальные линии в точках разрыва 
delta = 2    % задаем отступ от границ подынтерывалов
% в цикле находим значения полиномов и выводим их графики
for i = 1:8 
    X = ksi(i)-delta:0.01:ksi(i+1)+delta; % задаем значения аргумента
    P = polyval(C(i, :), X-ksi(i));  % вычисляем значения соответствующего полинома
    h = plot(X, P, 'LineWidth', 2);   % строим график полинома
    % под подынтервалом проводим горизонтальную линию того же цвета,
    % что и график полинома
    plot([ksi(i) ksi(i+1)], [-2 -2], 'Color', get(h, 'Color'), 'LineWidth', 6)  
end

Рис. 2. Полиномы, образующие сплайн.

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

% находим первую, вторую и третью производные сплайна
pp1 = fnder(pp);  
pp2 = fnder(pp1); % или pp2 = fnder(pp,2);
pp3 = fnder(pp2); % или pp3 = fnder(pp,3);
figure
% выводим график первой производной и вертикальные линии в точках разрыва
subplot(3, 1, 1)
fnplt(pp1)
hold on
plot([ksi; ksi], [-1 1], 'k:')
% выводим график второй производной и вертикальные линии в точках разрыва
subplot(3, 1, 2)
fnplt(pp2)
hold on
plot([ksi; ksi], [-1 1], 'k:')
% выводим график третьей производной и вертикальные линии в точках разрыва
subplot(3, 1, 3)
fnplt(pp3)
hold on
plot([ksi; ksi], [-1 1], 'k:')

Рис. 3. Первая, вторая и третья производные сплайна

Непрерывность первой и второй производных сплайна , а также разрывность его третьей производной можно установить также при помощи функции fnjmp, которая предназначена для нахождения скачков сплайна, например:

jms3 = fnjmp(pp3, ksi)

jms3 =
         0    0.4792   -1.6771    0.0369    1.5296    0.0369   -1.6771    0.4792         0

Нулевые скачки в первой и последней точках разрыва объясняются тем, что вне базового интервала сплайн является соответствующим полиномом: первым (для ) и последним (для ).

Итак, кубический сплайн и его производные до второго порядка включительно непрерывны во внутренних точках разрыва . Это накладывает на сплайн 7×3 = 21 условие. Еще 9 условий берутся из совпадения сплайна и функции cosx в точках . Итого, 21+ 9 = 30 условий. Однако, сплайн образован из 8 полиномов четвертого порядка (т.е. третьей степени). Следовательно, всего надо определить 8×4 = 32 коэффициента. Для единственности сплайна накладываются дополнительные условия в точках , называемые граничными условиями.

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

В нашем примере совпадает со значением производной кубического полинома, проходящего через точки , а совпадает со значением производной кубического полинома, проходящего через точки . В этом несложно убедиться на примере левой границы (см. рис. 4). Для построения полинома, проходящего через некоторые точки можно воспользоваться, например, функцией polyfit, входящей в MATLAB. Функция polyfit возвращает коэффициенты полинома заданной степени, приближающего данные в смысле наименьших квадратов. Однако, если число точек ровно на единицу превосходит степень, то такой полином будет как раз интерполяционным. Для вычисления производной полинома применим функцию MATLAB polyder, которая возвращает коэффициенты полинома, являющегося производной исходного.

% строим интерполяционный полином, проходящий через точки 
% (ksi(1), y(1)),   (ksi(2), y(2)),    (ksi(3), y(3)),    (ksi(4), y(4))
p3 = polyfit(ksi(1:4), y(1:4), 3)
% вычисляем значение интерполяционного полинома 
% в точках отрезка [ksi(1)-1.2, ksi(5)]
xi = ksi(1)-1.2:0.01:ksi(5);
yi = polyval(p3, xi);
% строим графики сплайна и полинома 
figure 
fnplt(pp, [ksi(1) ksi(5)])  % выводим график сплайна на отрезке [ksi(1) ksi(5)]
hold on
plot(xi, yi, 'g') % выводим график полинома на отрезке [ksi(1)-1.2, ksi(5)]
plot(ksi(1:4), y(1:4), 'ok') % отображаем маркерами точки разрыва 

% находим значение производной полинома в ksi(1)
p3dx = polyder(p3)
p3slope = polyval(p3dx, ksi(1));
% находим значение производной сплайна в ksi(1)
pp1 = fnder(pp); 
ppslope = fnval(pp1, ksi(1));
% выводим разность значений производных полинома и сплайна в ksi(1)
delta=p3slope-ppslope
% определяем анонимную функцию для касательной к сплайну в точке ksi(1)
ytangent = @(t) p3slope*(t-x(1))+y(1);
% строим график касательной
fplot(ytangent, [-1 3], 'k')
% устанавливаем пределы осей и размещаем легенду
set(gca, 'XLim', [ksi(1)-1.2 ksi(4)+0.2], 'YLim', [-1.2 1.5])
legend('spline', 'polynom', 'sites', 'tangent')

Рис. 4. Граничные условия, используемые csape по умолчанию

Примечание

Для конструирования интерполяционного полинома вовсе не обязательно было использовать функцию polyfit. Можно было снова обратиться к csape:

P3 = csape(ksi(1:4), y(1:4))

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

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


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

Система Orphus

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