В этом примере показано, как использовать команды от Curve Fitting Toolbox™, чтобы соответствовать сплайну к тестовым данным титана с ручным и автоматическим выбором узлов.
Вот некоторые данные, которые записывают определенное свойство титана, измеренного в зависимости от температуры. Мы будем использовать его, чтобы проиллюстрировать некоторые проблемы с интерполяцией сплайна.
[xx,yy] = titanium;
График данных показывает довольно резкий пик.
plot(xx,yy,'bx');
frame = [-10 10 -.1 .3]+[min(xx),max(xx),min(yy),max(yy)];
axis(frame);
Мы выбираем несколько точек данных из этих несколько грубых данных, поскольку мы хотим интерполировать. Вот изображение данных с выбранными отмеченными точками данных.
pick = [1 5 11 21 27 29 31 33 35 40 45 49]; tau = xx(pick); y = yy(pick); hold on plot(tau,y,'ro'); hold off
Поскольку сплайн порядка k с узлами n+k имеет n степени свободы, и у нас есть 12 точек данных, подгонка с четвертым сплайном порядка требует 12+4 = 16 узлов. Кроме того, эта последовательность узла t должна быть такова, что i-ый сайт данных находится в поддержку i-ого B-сплайна. Мы достигаем этого при помощи сайтов данных как узлы, но добавляем два простых узла с обоих концов.
dl = tau(2) - tau(1); dr = tau(end) - tau(end-1); t = [tau(1)-dl*[2 1] tau tau(end)+dr*[1 2]]; % construct the knot sequence plot(tau,y,'ro'); hold on axis(frame+[-2*dl 2*dr 0 0]) plot(t,repmat(frame(3)+.03,size(t)),'kx') hold off legend({'Data Values' 'Knots'},'location','NW')
Мы используем эту последовательность узла, чтобы создать интерполирующий кубический сплайн.
sp = spapi(t,tau,y);
Теперь для графика. Поскольку мы не заботимся о части сплайна вне интервала данных, мы ограничиваем график тем интервалом.
plot(tau,y,'ro') axis(frame) hold on fnplt(sp,[tau(1) tau(end)], 'k') hold off
Более внимательное рассмотрение в левой части подгонки сплайна показывает некоторые волнистости.
xxx = linspace(tau(1),tau(5),41); plot(xxx, fnval(sp, xxx), 'k', tau, y, 'ro'); axis([tau(1) tau(5) 0.6 1.2]);
Неблагоразумный удар в первом интервале происходит от того, что наш сплайн идет гладко, чтобы обнулить на уровне его первого узла. Видеть, что, вот изображение целого сплайна, наряду с его последовательностью узла и точками данных.
fnplt(sp,'k'); hold on plot(tau,y,'ro', t,repmat(.1,size(t)),'kx'); hold off legend({'Spline Interpolant' 'Data Values' 'Knots'},'location','NW')
Вот простой способ осуществить более разумное поведение на границе. Мы добавляем еще две точки данных вне интервала определенных данных и выбираем в качестве наших данных там значения прямой линии через первые две точки данных.
tt = [tau(1)-[4 3 2 1]*dl tau tau(end)+[1 2 3 4]*dr]; xx = [tau(1)-[2 1]*dl tau tau(end)+[1 2]*dr]; yy = [y(1)-[2 1]*(y(2)-y(1)) y y(end)+[1 2]*(y(end)-y(end-1))]; sp2 = spapi(tt,xx,yy); plot(tau,y,'ro', xx([1 2 end-1 end]),yy([1 2 end-1 end]),'bo'); axis(frame+[-2*dl 2*dr 0 0]); hold on fnplt(sp2,'b',tau([1 end])) hold off legend({'Original Data' 'Data Added for End Conditions' ... 'Fit with Added Data'},'location','NW')
Вот сравнение двух подгонок сплайна, чтобы показать сокращение волнистости в первом и последнем интервале.
hold on fnplt(sp,'k',tau([1 end])) hold off legend({'Original Data' 'Data Added for End Conditions' ... 'Fit with Added Data' 'Original Fit'},'location','NW')
Наконец, вот более внимательное рассмотрение в первых четырех интервалах данных, которое показывает более ясно сокращение волнистости около левого конца.
plot(tau,y,'ro',xxx,fnval(sp2,xxx),'b',xxx,fnval(sp,xxx),'k'); axis([tau(1) tau(5) .6 1.2]); legend({'Original Data' 'Fit with Added Data' ... 'Original Fit'},'location','NW')
Если вся эта деталь выключает вас, позвольте Curve Fitting Toolbox выбрать узлы для вас. Задайте желаемый порядок interpolant как первый входной параметр к команде интерполяции сплайна spapi, вместо последовательности узла.
autosp = spapi(4, tau, y); knots = fnbrk(autosp,'knots'); plot(tau, y, 'ro') hold on fnplt(autosp,'g') plot(knots, repmat(.5,size(knots)),'gx') hold off legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ... 'Knots Chosen by SPAPI'}, 'location','NW')
Ниже результат намного лучшего выбора узла, полученного путем сдвига узла в 842 немного направо и узла в 985 немного налево.
knots([7 12]) = [851, 971]; adjsp = spapi(knots, tau, y); hold on fnplt(adjsp,'r',2) plot(knots, repmat(.54,size(knots)),'rx') hold off legend({'Data Values' 'Fit With Knots Chosen by SPAPI' ... 'Knots Chosen by SPAPI' 'Fit With Knots Adjusted' ... 'Adjusted Knots'}, 'location','NW')
Или иначе, просто попробуйте стандартный кубический сплайн interpolant, предоставленный csapi. Это составляет немного отличающийся выбор узлов.
autocs = csapi(tau, y); plot(tau, y, 'ro') hold on fnplt(autocs,'c') hold off
С такими быстро различными данными трудно получить соглашение среди всего разумного interpolants, даже если каждый из них является кубическим сплайном. График ниже показов все пять interpolants, для сравнения.
plot(tau, y, 'ro') hold on fnplt(sp,'k',tau([1 end])) % black: original fnplt(sp2,'b',tau([1 end])) % blue: with special end conditions fnplt(autosp,'g') % green: automatic knot choice by SPAPI fnplt(autocs,'c') % cyan: automatic knot choice by CSAPI fnplt(adjsp,'r',2) % red: knot choice by SPAPI slightly changed hold off legend({'Data Values' 'Original Fit' 'Special End Conditions' ... 'With Knots Chosen by SPAPI' 'With Knots Chosen by CSAPI' ... 'With Adjusted Knots'},'location','NW')
