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

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

Шлицуйте Interpolants различных порядков сглаженных данных

Даже кубический сплайн 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

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

Вторая производная двух кубических сплайнов Interpolants тех же сглаженных данных

Поскольку вторая производная кубического сплайна является прерывистой линией с вершинами в пропусках сплайна, вы видите ясно, что 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);
Для просмотра документации необходимо авторизоваться на сайте