Как выбрать узлы

Этот пример показывает, как выбрать и оптимизировать узлы с помощью 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

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

Используя 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

Это довольно хорошо. Кстати, даже один внутренний узел меньше не были бы достаточны в этом случае.

Для просмотра документации необходимо авторизоваться на сайте