Документация

Кубический сплайн Interpolant сглаженных данных

Предположим, что вы хотите интерполировать некоторые сглаженные данные, например, к

rng(6), x = (4*pi)*[0 1 rand(1,15)]; y = sin(x);

Можно использовать кубический сплайн interpolant полученный

cs = csapi(x,y);

и постройте сплайн, наряду с данными, со следующим кодом:

fnplt(cs);
hold on
plot(x,y,'o')
legend('cubic spline','data')
hold off

Это производит фигуру как следующее.

Это - более точно, кубический сплайн interpolant с граничными условиями не-узла, означая, что это - уникальный кусочный кубический полином с двумя непрерывными производными с пропусками на всех внутренних сайтах данных за исключением крайнего левого и самого правого. Это - тот же interpolant, как произведено командой MATLAB^® spline, spline(x,y).

Периодические данные

Синусоидальная функция 2π-periodic. Чтобы проверять, как хорошо ваш interpolant делает в этом отношении, вычислите, например, различие в значении его первой производной в этих двух конечных точках,

diff(fnval(fnder(cs),[0 4*pi]))
ans = -.0100

который не так хорош. Если вы предпочитаете получать interpolant, первые и вторые производные которого в этих двух конечных точках, 0 и 4*pi, соответствии, используют вместо этого команду csape, который разрешает спецификацию многих различных видов граничных условий, включая периодические граничные условия. Так, используйте вместо этого

pcs = csape(x,y,'periodic');

для которого вы добираетесь

diff(fnval(fnder(pcs),[0 4*pi]))

Выводом является ans = 0 как различие наклонов конца. Даже различие в конце вторые производные является небольшим:

diff(fnval(fnder(pcs,2),[0 4*pi]))

Выводом является ans = -4.6074e-015.

Другие граничные условия

Другие граничные условия могут быть обработаны также. Например,

cs = csape(x,[3,y,-4],[1 2]);

предоставляет кубическому сплайну interpolant пропуски в и его наклон на крайнем левом сайте данных, равном 3, и его вторая производная на самом правом сайте данных, равном-4.

Общая интерполяция сплайна

Если вы хотите интерполировать на сайтах кроме пропусков и/или сплайнами кроме кубических сплайнов с простыми узлами, то вы используете команду spapi. В его самой простой форме вы сказали бы sp = spapi(k,x,y); в котором первый аргумент, k, задает порядок сплайна интерполяции; это - количество коэффициентов в каждой полиномиальной части, т.е. еще 1, чем номинальная степень ее полиномиальных частей. Например, следующие данные показывают линейное, квадратичное, и биквадратный сплайн interpolant к вашим данным, как получено операторами

sp2 = spapi(2,x,y); fnplt(sp2,2), hold on
sp3 = spapi(3,x,y); fnplt(sp3,2,'k--'), 
sp5 = spapi(5,x,y); fnplt(sp5,2,'r-.'), plot(x,y,'o')
legend('linear','quadratic','quartic','data'), hold off

Даже кубический сплайн interpolant полученный из spapi отличается от того, обеспеченного csapi и spline. Чтобы подчеркнуть их различие, вычислите и постройте их вторые производные, можно следующим образом:

fnplt(fnder(spapi(4,x,y),2)), hold on,
fnplt(fnder(csapi(x,y),2),2,'k--'),plot(x,zeros(size(x)),'o')
legend('from spapi','from csapi','data sites'), hold off

Это дает следующий график:

Поскольку вторая производная кубического сплайна является прерывистой линией с вершинами в пропусках сплайна, вы видите ясно, что csapi помещает пропуски на сайтах данных, в то время как spapi не делает.

Свяжите выбор узлом

На самом деле, возможно задать явным образом, где сплайн interpolant должен иметь свои пропуски, с помощью команды sp = spapi(knots,x,y); в котором последовательность предоставления knots, определенным способом, пропуски, которые будут использоваться. Например, вспоминание, что вы выбрали y, чтобы быть sin(x), командой

ch = spapi(augknt(x,4,2),[x x],[y cos(x)]);

предоставляет кубического Эрмита interpolant синусоидальной функции, а именно, кусочная кубическая функция, с пропусками во всем x(i), который совпадает с синусоидальной функцией в значении и наклоне во всем x(i). Это делает interpolant непрерывное с непрерывной первой производной, но в целом она имеет скачки через перерывы в ее второй производной. Как эта команда знает, какая часть массива значения данных  [y cos(x)] предоставляет значения и который наклоны? Заметьте, что массив сайта данных здесь дан как [x x], т.е. каждый сайт данных появляется дважды. Также заметьте, что y(i) сопоставлен с первым вхождением x(i), и cos(x(i)) сопоставлен со вторым вхождением x(i). Значение данных, сопоставленное с первым выступлением сайта данных, взято, чтобы быть значением функции; значение данных, сопоставленное со вторым внешним видом, взято, чтобы быть наклоном. Если бы был третий внешний вид того сайта данных, соответствующее значение данных было бы взято в качестве второго производного значения, которое будет соответствующим на том сайте. Смотрите Построение и Работу со Сплайнами B-формы для обсуждения команды, augknt раньше здесь генерировал соответствующую "последовательность узла".

Сглаживание

Что, если данные являются шумными? Например, предположите, что данные значения

noisy = y + .3*(rand(size(x))-.5);

Затем вы можете предпочесть аппроксимировать вместо этого. Например, вы можете попробовать кубический сплайн сглаживания, полученный командой

scs = csaps(x,noisy);

и построенный по

fnplt(scs,2), hold on, plot(x,noisy,'o'),
legend('smoothing spline','noisy data'), hold off

Это производит фигуру как это:

Если вам не нравится уровень сглаживания сделанного csaps(x,y), можно изменить его путем определения параметра сглаживания, p, в качестве дополнительного третьего аргумента. Выберите этот номер где угодно между 0 и 1. Когда p изменяется с 0 до 1, изменения сплайна сглаживания, соответственно, от одного экстремального значения, наименьшие квадраты прямолинейное приближение к данным, к другому экстремальному значению, "естественный" кубический сплайн interpolant к данным. Поскольку csaps возвращает параметр сглаживания, на самом деле используемый в качестве дополнительного второго вывода, вы могли теперь экспериментировать, можно следующим образом:

[scs,p] = csaps(x,noisy); fnplt(scs,2), hold on
fnplt(csaps(x,noisy,p/2),2,'k--'), 
fnplt(csaps(x,noisy,(1+p)/2),2,'r:'), plot(x,noisy,'o')
legend('smoothing spline','more smoothed','less smoothed',...
'noisy data'), hold off

Это производит следующее изображение.

Время от времени вы можете предпочесть просто получать самый сглаженный кубический сплайн sp, который является в заданном допуске tol определенных данных в том смысле, что norm(noisy - fnval(sp,x))^2 <= tol. Вы создаете этот сплайн с командой sp = spaps(x,noisy,tol) для вашего заданного допуска tol.

Наименьшие квадраты

Если вы предпочитаете аппроксимирующую функцию наименьших квадратов, можно получить ее оператором sp = spap2(knots,k,x,y); в котором и последовательность узла knots и порядок должен быть обеспечен k сплайна.

Популярный выбор для порядка равняется 4, и это дает вам кубический сплайн. Если у вас нет четкого представления о том, как выбрать узлы, просто задать количество полиномиальных частей, вы хотите используемый. Например,

sp = spap2(3,4,x,y); 

дает кубический сплайн, состоящий из трех полиномиальных частей. Если получившаяся ошибка неровна, вы можете попробовать за лучшее распределение узла при помощи newknt можно следующим образом:

sp = spap2(newknt(sp),4,x,y);