Конструкция и работа с PPFORM

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

Введение

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

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

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

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

Вот пример pp, состоящего из трех квадратичных полиномов, т.е. 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

Точное описание pp в терминах последовательности пропуска 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.

Для pp на рисунке выше, 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) определяется расширением первой или последней полиномиальной части.

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

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, полное описание этой функции pp в так называемом ppform.

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

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

Чтобы вычислить pp, используйте 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

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

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

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

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

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, вторая производная там не бесконечна. Это имеет последствия, когда мы интегрируем вторую производную.

Чтобы интегрировать pp, используйте 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, который не может быть восстановлен, поскольку интеграл любой функции pp непрерывен.

Получить детали можно с помощью команды fnbrk. Для примера

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

     0     1     4     6

восстанавливает пропуск последовательность pp в 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])

Векторные Кусочные полиномы

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

Для примера вот pp с 2 векторов значениями, описывающий модуль квадрат как его показы график. Это 2D-curve, поэтому его размерность равна 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')

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

pp в Curve Fitting Toolbox также может быть многомерным, а именно, тензорным продуктом одномерных функций pp. ppform такой многомерной pp лишь немного сложнее, с breaks теперь массив ячеек, содержащий последовательность пропусков для каждой переменной, и coefs теперь многомерный массив. Намного сложнее составить неслучайную такую функцию с нуля, поэтому здесь мы не будем стараться, тем более, что тулбокс предназначен для помощи в конструкции функций pp из некоторых условий о них. Например, сфера на этом рисунке построена с помощью 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-форму, то есть его представление как линейной комбинации B-сплайнов. Для получения дополнительной информации смотрите Конструкция и Работа с B-формой.