Чебышев (a.k.a. Equioscillating) заданный сплайн

По определению, для данной последовательности узла t из длины n+k, C = C_{t,k} уникальный элемент S_{t,k} из макс. нормы 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');

Итерация

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

Нахождение их экстремальное значение является самостоятельно итеративным процессом, а именно, для нахождения нулей производного 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);

Новое приближение является более близко сплайном equioscillating.

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));

как свидетельствуется близостью equioscillation ошибки.

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