Построение и работа с PPFORM

В этом примере показано, как создать и работать с ppform сплайна в Curve Fitting Toolbox™.

Введение

(Одномерный) кусочный полином или стр для краткости характеризуется его последовательностью пропуска, breaks скажите, и его массив коэффициентов, coefs скажите локальной формы степени ее полиномиальных частей. Последовательность пропуска принята, чтобы строго увеличиться,

  breaks(1) < breaks(2) < ... < breaks(l+1),

с l количество полиномиальных частей, которые составляют стр. В рисунке ниже, breaks [0,1,4,6], следовательно l 3.

В то время как эти полиномы могут иметь различные степени, они все зарегистрированы как полиномы того же порядка k, т.е. массив коэффициентов coefs имеет размер [l,k], с coefs(j,:) содержа k коэффициенты в локальной степени формируются для j- часть полинома th.

Вот пример, стр составили из трех квадратичных полиномов, т.е. l = k = 3. Пропуски отмечены красным.

sp = spmak([0 1 4 4 6],[2 -1]);
pp = fn2fm(sp,'pp') ;
breaks = fnbrk(pp,'b');
coefs = fnbrk(pp,'c');
coefs(3,[1 2]) = [0 1];
pp = ppmak(breaks,coefs,1);
fnplt(pp,[breaks(1)-1 breaks(2)],'g',1.8)
hold on
fnplt(pp, breaks([2 3]),'b',1.8)
fnplt(pp,[breaks(3),breaks(4)+1],'m',1.8)
lp1 = length(breaks);
xb = repmat(breaks,3,1);
yb = repmat([2;-2.2;NaN],1,lp1);
plot(xb(:),yb(:),'r')
axis([-1 7 -2.5 2.3])
hold off

Точное описание стр в терминах последовательности пропуска breaks и массив коэффициентов coefs

  pp(t) = polyval(coefs(j,:), t-breaks(j))
                              for breaks(j) <= t < breaks(j+1)

где, чтобы вспомнить,

  polyval(a,x) = a(1)*x^(k-1) + a(2)*x^(k-2) + ... + a(k)*x^0.

Для стр на рисунке выше, breaks(1) 0, и coefs(1,:) [-1/2 0 0], в то время как breaks(3) 4, и coefs(3,:) [0 1 - 1]. Для t не в [breaks(1) .. breaks(l+1)], pp(t) задан путем расширения первой или последней полиномиальной части.

Стр обычно создается посредством процесса интерполяции или приближения. Но также возможно составить один в ppform с нуля, с помощью команды ppmak. Например, стр выше могут быть получены как

breaks = [0 1 4 6];
coefs = [1/2 0 0 -1/2 1 1/2 0 1 -1];
fn = ppmak(breaks,coefs)
fn = 

  struct with fields:

      form: 'pp'
    breaks: [0 1 4 6]
     coefs: [3x3 double]
    pieces: 3
     order: 3
       dim: 1

Это возвращается в fn, полное описание этого стр функционирует в так называемой ppform.

Различные команды в Curve Fitting Toolbox могут работать с этой формой. Остающиеся разделы показывают некоторые примеры.

Работа на кусочных полиномах

Чтобы оценить стр, используйте fnval команда.

fnval(fn, -1:7)
ans =

  Columns 1 through 7

    0.5000         0    0.5000    1.0000    0.5000   -1.0000         0

  Columns 8 through 9

    1.0000    2.0000

Чтобы дифференцировать стр, используйте fnder команда.

dfn = fnder ( fn );
hold on
fnplt(dfn, 'jumps','y', 3)
hold off
h1 = findobj(gca,'Color','y');
legend(h1,{'First Derivative'},'location','SW')

Обратите внимание на то, что производная стр в качестве примера непрерывна в 1, но прерывиста в 4. Также обратите внимание что, по умолчанию, fnplt строит ppform на ее основном интервале, т.е. на интервале [breaks(1) .. breaks(end)].

Можно также использовать fnder взять вторую производную стр.

ddfn = fnder(fn, 2);
hold on
fnplt( ddfn ,'j', 'k', 1.6)
hold off
h2 = findobj(gca,'Color','k');
legend([h1 h2],{'First Derivative' 'Second Derivative'},'location','SW')

Вторая производная является кусочной константой.

Обратите внимание на то, что дифференцирование через fnder сделан отдельно для каждой полиномиальной части. Например, несмотря на то, что первая производная имеет разрыв скачка через 4, вторая производная весьма конечна там. Это имеет последствия, когда мы интегрируем вторую производную.

Чтобы интегрировать стр, используйте fnint команда.

iddfn = fnint(ddfn);
hold on
fnplt(iddfn, 'b', .5)
hold off
h3 = findobj(gca,'Color','b', 'LineWidth',.5);
legend([h1 h2 h3],{'First Derivative' 'Second Derivative' ...
                   'Integral of Second Derivative'},'location','SW')

Обратите внимание на то, что интегрирование второй производной действительно восстанавливает первую производную, за исключением скачка через 4, который не может быть восстановлен, поскольку интеграл любой функции стр непрерывен.

Можно получить части при помощи команды fnbrk. Например,

breaks = fnbrk(fn, 'breaks')
breaks =

     0     1     4     6

восстанавливает последовательность пропуска стр в fnв то время как

piece2 = fnbrk(fn, 2);

восстанавливает вторую полиномиальную часть, как этот график подтверждает.

fnplt(pp,[breaks(1)-1 breaks(2)],'g',1.8)
hold on
fnplt(piece2, 'b', 2.5, breaks([2 3])+[-1 .5])
fnplt(pp,[breaks(3),breaks(4)+1],'m',1.8)
plot(xb(:),yb(:),'r')
title('The Polynomial that Supplies the Second Polynomial Piece')
hold off
axis([-1 7 -2.5 2.3])

Кусочные полиномы с векторным знаком

Стр могут также быть с векторным знаком, чтобы описать кривую, в с 2 пробелами или с 3 пробелами. В этом случае каждый локальный полиномиальный коэффициент является вектором, а не номером, но ничем иным об изменениях ppform. Существует одна дополнительная часть ppform, чтобы записать это, размерность ее цели.

Например, вот оцененные стр 2 векторов, описывающие модульный квадрат, когда его график показывает. Это - 2D кривая, следовательно ее размерность равняется 2.

square = ppmak(0:4, [1 0  0 1  -1 1  0 0 ; 0 0  1 0  0 1  -1 1]);
fnplt(square,'r',2)
axis([-.5 1.5 -.5 1.5])
axis equal
title('A Vector-Valued PP that Describes a Square')

Многомерные кусочные полиномы

Стр в Curve Fitting Toolbox могут также быть многомерными, а именно, продукт тензора одномерных функций стр. ppform таких многомерных стр незначительно более сложна с breaks теперь массив ячеек, содержащий последовательность пропуска для каждой переменной и coefs теперь многомерный массив. Намного более трудно составить неслучайное такая функция с нуля, таким образом, мы не попробуем это здесь, особенно поскольку тулбокс предназначается, чтобы помочь с конструкцией функций стр от некоторых условий о них. Например, сфера в этом рисунке создается при помощи csape.

x = 0:4;
y = -2:2;
s2 = 1/sqrt(2);
v = zeros(3,7,5);
v(3,:,:) = [0 1 s2 0 -s2 -1 0].'*[1 1 1 1 1];
v(2,:,:) = [1 0 s2 1 s2 0 -1].'*[0 1 0 -1 0];
v(1,:,:) = [1 0 s2 1 s2 0 -1].'*[1 0 -1 0 1];
sph = csape({x,y},v,{'clamped','periodic'});
fnplt(sph)
axis equal
axis off
title('A Sphere Described by a Bicubic 3-Vector-Valued Spline')

В то время как ppform сплайна эффективна для оценки, конструкция сплайна из некоторых данных обычно обрабатывается более эффективно первым определением его B-формы, i.e., его представление как линейная комбинация B-сплайнов. Для получения дополнительной информации смотрите Построение и работу с B-формой.