Этот пример показывает, как выбрать и оптимизировать узлы с помощью optknt
и команд newknt
от Curve Fitting Toolbox™.
Вот некоторые выборочные данные, очень используемые для тестирования приближения сплайна с переменными узлами, так называемыми Данными о Тепле Титана. Они записывают некоторое свойство титана, измеренного как функция температуры.
[xx,yy] = titanium; plot(xx,yy,'x'); axis([500 1100 .55 2.25]); title('The Titanium Heat Data'); hold on
Заметьте довольно резкий пик. Мы будем использовать эти данные, чтобы проиллюстрировать некоторые методы для выбора узла.
Во-первых, мы выбираем несколько точек данных от этих несколько грубых данных. Мы интерполируем использование этого подмножества, затем сравним результаты с полным набором данных.
pick = [1 5 11 21 27 29 31 33 35 40 45 49]; tau = xx(pick); y = yy(pick); plot(tau,y,'ro'); legend({'Full Dataset' 'Subsampled Data'}, 'location','NW');
Сплайн порядка k
с узлами n+k
имеет степени свободы n
. Поскольку у нас есть 12 сайтов данных, tau(1) < ... < tau(12)
, подгонка с кубическим сплайном, т.е. четвертым сплайном порядка, требует последовательности узла t
длины 12+4.
Кроме того, последовательность узла t
должен удовлетворить условия Шенберга-Уитни, т.е. должен быть таков, что i-th сайт данных находится в поддержку i-th B-сплайна, т.е.
t(i) < tau(i) < t(i+k) for all i,
с равенством, позволенным только в случае узла кратности k
.
Один способ выбрать последовательность узла, удовлетворяющую все эти условия, как оптимальные узлы Gaffney/Powell и Micchelli/Rivlin/Winograd.
В оптимальной интерполяции сплайна, к значениям на сайтах
tau(1), ..., tau(n)
скажите, узлы выбраны, чтобы минимизировать константу в формуле стандартной погрешности. А именно, первое и последний сайт данных выбраны в качестве узлов k-сгиба. Остающиеся узлы n-k
предоставляются optknt
.
Вот начало справки от optknt
:
OPTKNT Оптимальное распределение узла.
OPTKNT(TAU,K) returns an `optimal' knot sequence for
interpolation at data sites TAU(1), ..., TAU(n) by splines of
order K. TAU must be an increasing sequence, but this is not
checked.
OPTKNT(TAU,K,MAXITER) specifies the number MAXITER of iterations
to be tried, the default being 10.
The interior knots of this knot sequence are the n-K
sign-changes in any absolutely constant function h ~= 0 that
satisfies
integral{ f(x)h(x) : TAU(1) < x < TAU(n) } = 0
for all splines f of order K with knot sequence TAU.
Мы пытаемся использовать optknt
для интерполяции на нашем примере, интерполируя кубическими сплайнами к данным
(tau(i), y(i)), for i = 1, ..., n.
k = 4; osp = spapi( optknt(tau,k), tau,y); fnplt(osp,'r'); hl = legend({'Full Dataset' 'Subsampled Data' ... 'Cubic Spline Interpolant Using Optimal knots'}, ... 'location','NW'); hl.Position = hl.Position-[.14,0,0,0];
Это является немного дезориентирующим!
Здесь, отмеченный звездами, также внутренние оптимальные узлы:
xi = fnbrk(osp,'knots'); xi([1:k end+1-(1:k)]) = []; plot(xi,repmat(1.4, size(xi)),'*'); hl = legend({'Full Dataset' 'Subsampled Data' ... 'Cubic Spline Interpolant Using Optimal knots' ... 'Optimal Knots'}, 'location','NW'); hl.Position = hl.Position-[.14,0,0,0];
Выбор узла для оптимальной интерполяции разработан, чтобы сделать максимум по всем функциям f
отношения
norm(f - If) / norm(D^k f)
как можно меньше, где числитель является нормой ошибки интерполяции, f - If
, и знаменатель является нормой k
-th производная interpolant, D^k f
. Поскольку наши данные подразумевают, что D^k f
является довольно большим, ошибка интерполяции около плоской части данных имеет приемлемый размер для такой 'оптимальной' схемы.
На самом деле, для этих данных, обычный кубический сплайн interpolant обеспеченный csapi
вполне успевает:
cs = csapi(tau,y); fnplt(cs,'g',2); hl = legend({'Full Dataset' 'Subsampled Data' ... 'Cubic Spline Interpolant Using Optimal knots' ... 'Optimal Knots' 'Cubic Spline Interpolant Using CSAPI'}, ... 'location','NW'); hl.Position = hl.Position-[.14,0,0,0]; hold off
Узлы должны быть выбраны при выполнении приближения наименьших квадратов сплайнами. Один подход должен использовать равномерно распределенные узлы для начала, затем использовать newknt
с приближением, полученным для лучшего распределения узла.
Следующие разделы иллюстрируют, что эти шаги с полным титаном нагревают набор данных.
Мы запускаем с универсальной последовательности узла.
unif = linspace(xx(1), xx(end), 2+fix(length(xx)/4)); sp = spap2(augknt(unif, k), k, xx, yy); plot(xx,yy,'x'); hold on fnplt(sp,'r'); axis([500 1100 .55 2.25]); title('The Titanium Heat Data'); hl = legend({'Full Dataset' ... 'Least Squares Cubic Spline Using Uniform Knots'}, ... 'location','NW'); hl.Position = hl.Position-[.14,0,0,0];
Это нисколько не удовлетворительно. Таким образом, мы используем newknt
для приближения сплайна того же порядка и с тем же количеством полиномиальных частей, но пропусками, лучше распределенными.
spgood = spap2(newknt(sp), k, xx,yy); fnplt(spgood,'g',1.5); hl = legend({'Full Dataset' ... 'Least Squares Cubic Spline Using Uniform Knots' ... 'Least Squares Cubic Spline Using NEWKNT'}, ... 'location','NW'); hl.Position = hl.Position-[.14,0,0,0]; hold off
Это довольно хорошо. Кстати, даже один внутренний узел меньше не были бы достаточны в этом случае.