exponenta event banner

Конструирование и работа с PPFORM

В этом примере показано, как создавать и работать с ppform сплайна в 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.

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

Работа с кусочными многочленами

Чтобы оценить 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 для записи этого, измерение его цели.

Например, вот 2-векторное значение pp, описывающее единичный квадрат, как показывает его график. Это 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 в панели инструментов фитинга кривой также может быть многомерным, а именно тензорным произведением одномерных 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-формой.