Чебышёв (а.к.а. Эквиоскилляция) Сплайн Определен

По определению для заданной последовательности узлов 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, имеет максимальную норму 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);

Затем мы используем два шага метода secant, получая итераты 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.'});