Чебышёв (а.к.а. Equioscillating) Определение сплайна

По определению, для данной узловой последовательности t длины n+k, C = C_{t,k} является уникальным элементом S_{t,k} из max-norm 1, который максимально колеблется на интервале [t_k .. t_{n+1}] и является положительным вблизи t_{n+1}. Это означает, что существует уникальный строго увеличивающийся tau длины n чтобы функция C в S_{k,t} дано

C(tau(i)) = (-1)^{n-i},

для всех i, имеет max-norm 1 на [t_k .. t_{n+1}]. Это подразумевает, что

tau(1) = t_k,

tau(n) = t_{n+1},

и что

t_i < tau(i) < t_{k+i},

для всех i. Действительно,

t_{i+1} <= tau(i) <= t_{i+k-1},

для всех i. Это приводит к тому, что последовательность узлов t предполагается, что такое неравенство возможно, что оказывается эквивалентным наличию всех элементов S_{k,t} непрерывный.

t = augknt([0 1 1.1 3 5 5.5 7 7.1 7.2 8], 4 );
[tau,C] = chbpnt(t,4);
xx = sort([linspace(0,8,201),tau]);
plot(xx,fnval(C,xx),'LineWidth',2);
hold on
breaks = knt2brk(t); 
bbb = repmat(breaks,3,1);
sss = repmat([1;-1;NaN],1,length(breaks));
plot(bbb(:), sss(:),'r');
hold off
ylim([-2 2]);
title('The Chebyshev Spline for a Particular Knot Sequence');
legend({'Chebyshev Spline' 'Knots'});

Короче говоря, шлиц Чебышёва C выглядит как многочлен Чебышёва. Выполняет аналогичные функции. Например, его крайность tau являются особенно хорошими сайтами для интерполяции из S_{k,t} так как норма результирующего проектора примерно так мала, как может быть.

hold on
plot(tau,zeros(size(tau)),'k+');
hold off
legend({'Chebyshev Spline' 'Knots' 'Extrema'});

Выбор сплайнового пространства

В этом примере мы пытаемся построить C для заданного сплайнового пространства.

Мы имеем дело с кубическими шлицами с простыми внутренними узлами, указанными

k = 4;
breaks = [0 1 1.1 3 5 5.5 7 7.1 7.2 8];
t = augknt(breaks, k)

t = 1×16

         0         0         0         0    1.0000    1.1000    3.0000    5.0000    5.5000    7.0000    7.1000    7.2000    8.0000    8.0000    8.0000    8.0000

таким образом получают сплайновое пространство размера

n = length(t)-k

n = 12

Начальное предположение

В качестве нашего первоначального предположения для tau, мы используем усредненные узловые значения

tau(i) = (t_{i+1} + ... + t_{i+k-1})/(k-1)

рекомендуется как хороший выбор сайта интерполяции и построить график полученного первого приближения к C.

tau = aveknt(t,k)

tau = 1×12

         0    0.3333    0.7000    1.7000    3.0333    4.5000    5.8333    6.5333    7.1000    7.4333    7.7333    8.0000

b = (-ones(1,n)).^(n-1:-1:0);
c = spapi(t,tau,b);
plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k');
hold on
fnplt(c,'r',1);
hold off
ylim([-2 2]);
title('First Approximation to an Equioscillating Spline');

Повторение

Для полного нивелирования используем алгоритм Ремеза. Это означает, что мы создаем новый tau как крайность нашего текущего приближения, cКому C и повторите попытку.

Нахождение этих крайностей само по себе является итеративным процессом, а именно, для нахождения нулей производной Dc нашего текущего приближения c.

Dc = fnder(c);

Берем нули управляющего многоугольника Dc как наше первое предположение для нулей Dc. Этот управляющий многоугольник имеет вершины (tstar(i),coefs(i)), где coefs являются коэффициентами Dc, и tstar узел усредняется.

[knots,coefs,np,kp] = fnbrk(Dc, 'knots', 'coefs', 'n', 'order');
tstar = aveknt(knots,kp);

Поскольку управляющий многоугольник является кусочно-линейным, его нули легко вычислить. Вот эти нули.

npp = 1:np-1;
guess = tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs));
fnplt(Dc,'r');
hold on
plot(tstar,coefs,'k.:');
plot(guess,zeros(1,np-1),'o');
hold off
title('First Derivative of the Approximation');
legend({'Dc' 'Control Polygon' 'Zeros of Control Polygon'});

Это дает очень хорошее первое предположение для фактических нулей Dc.

Теперь оцениваем Dc на обоих этих наборах участков.

sites = [guess; tau(2:n-1)];
values = fnval(Dc,sites);

Затем мы используем два шага секущего метода, получая итерации sites(3,:) и sites(4,:), с values(3,:) и values(4,:) соответствующие значения Dc.

sites(3:4,:) = 0;
values(3:4,:) = 0;
for j = 2:3
   rows = [j,j-1]; 
   Dcd = diff(values(rows,:));
   Dcd(Dcd==0) = 1; % guard against division by zero
   sites(j+1,:) = sites(j,:)-values(j,:).*(diff(sites(rows,:))./Dcd);
   values(j+1,:) = fnval(Dc,sites(j+1,:));
end

Мы принимаем последний итерат как вычисленные нули Dc, т.е. крайность нашего текущего приближения, c. Это наше новое предположение для tau.

tau = [tau(1) sites(4,:) tau(n)]

tau = 1×12

         0    0.2759    0.9082    1.7437    3.0779    4.5532    5.5823    6.5843    7.0809    7.3448    7.7899    8.0000

plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k');
hold on
fnplt(c,'r',1);
plot(guess,zeros(1,np-1),'o');
plot(tau(2:n-1),zeros(1,n-2),'x');
hold off
title('First Approximation to an Equioscillating Spline');
ax = gca;
h = ax.Children;
legend(h([3 1 2]),{'Approximation' 'Extrema' ...
                   'Zeros of First Derivative''s Control Polygon'});
axis([0 8 -2 2]);

Окончание первого шага итерации

Вычисляем полученное новое приближение к сплайну Чебышева, используя новое предположение для tau.

cnew = spapi(t,tau,b);

Новая аппроксимация является более близким к эквиосциллирующему сплайну.

plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k');
hold on
fnplt(c,'r',1);
fnplt(cnew, 'k', 1);
hold off
ax = gca;
h = ax.Children;
legend(h([2 1]),{'First Approximation' 'Updated Approximation'});
axis([0 8 -2 2]);

Если это недостаточно близко, просто повторите попытку, начиная с этого нового tau. Для этого конкретного примера следующая итерация уже предоставляет чебышёвский сплайн для графической точности.

Использование очков Чебышева-Демко

Сплайн Чебышева для заданного сплайнового пространства S_{k,t}, наряду с его экстремальными значениями, доступны в качестве необязательных выходов из chbpnt на панели инструментов. Эти крайности были предложены в качестве хороших интерполированных участков Стивеном Демко, отсюда теперь называются участками Чебышева-Демко. В этом разделе показан пример их использования.

Если вы решили аппроксимировать функцию квадратного корня на интервале [0 .. 1] кубическими шлицами с узловой последовательностью

k = 4;
n = 10;
t = augknt(((0:n)/n).^8,k);

тогда хорошее приближение к функции квадратного корня из этого определенного сплайнового пространства дается

tau = chbpnt(t,k);
sp = spapi(t,tau,sqrt(tau));

о чем свидетельствует близкое эквиосцилирование ошибки.

xx = linspace(0,1,301);
plot(xx, fnval(sp,xx)-sqrt(xx));
title({'Error in Interpolant to Square Root','at Chebyshev-Demko Sites.'});