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