Предположим, что вы хотите интерполировать некоторые сглаженные данные, e.g., к
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 делает в этом отношении, вычислите, e.g., различие в значении его первой производной в этих двух конечных точках,
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, задает порядок сплайна интерполяции; это - количество коэффициентов в каждой полиномиальной части, i.e., еще 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);