pp 型の作成と操作

はじめに

(一変量) 区分的多項式 (略して pp) は、その "ブレーク シーケンス" (breaks など) と多項式区分の局所的な指数形式の "係数配列" (coefs など) によって指定されます。ブレーク シーケンスは厳密に増加するものと仮定します。

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

l は、pp を構成する多項式区分の数です。次の Figure で、breaks は [0,1,4,6] であるため、l は 3 です。

これらの多項式の次数はさまざまですが、すべて同じ "次数" k の多項式として記録されます。つまり、係数配列 coefs はサイズ [l,k] で、coefs(j,:) には j 番目の多項式区分の局所的な指数形式で、最大から最小の指数の順に、k 係数が含まれます。

次に、3 つの 2 次多項式 (すなわち l = k = 3) で構成される PP の例を示します。ブレークは赤でマークされています。

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 の正確な記述は次のとおりです。

  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.

上の Figure の PP の場合、breaks(1) は 0 で coefs(1,:) は [-1/2 0 0] ですが、breaks(3) は 4 で coefs(3,:) は [0 1 -1] です。[breaks(1) .. breaks(l+1)] に含まれない t の場合、pp(t) は最初または最後の多項式区分を拡張することによって定義されます。

通常、pp は内挿または近似のプロセスを通じて作成されます。ただし、コマンド ppmak を使用してゼロから pp 型を構成することもできます。たとえば、上記の 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: [3×3 double]
    pieces: 3
     order: 3
       dim: 1

これは、いわゆる pp 型でこの pp 関数の詳細説明を fn で返します。

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 により "基本区間"、すなわち区間 [breaks(1) .. breaks(end)] 上に pp 型がプロットされることにも注意してください。

fnder を使用して、pp の 2 階微分を取得することもできます。

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

この 2 階微分は区分的定数となります。

fnder による微分は多項式区分ごとに別々に行うことに注意してください。たとえば、1 階微分には 4 にまたがるジャンプ不連続点がありますが、そこでの 2 階微分は無限ではありません。これは 2 階微分を積分する際に影響があります。

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

任意の pp 関数の積分は連続であるため、2 階微分を積分すると、復元されない 4 をまたがるジャンプを除いては 1 階微分が復元されることに注意してください。

コマンド fnbrk により、部分を取得できます。たとえば、

breaks = fnbrk(fn, 'breaks')

breaks =

     0     1     4     6

により、fn 内の pp のブレーク シーケンスが復元されますが、一方、

piece2 = fnbrk(fn, 2);

では、このプロットで確認されるように 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 次元空間の曲線を記述することもできます。その場合、局所的な各多項式係数は数値ではなくベクトルとなりますが、pp 型に関して他に変わるものはありません。これを記録するための pp 型の追加部分として、そのターゲットの "次元" があります。

例として、次に単位正方形を表す 2 ベクトル値の pp を示します (プロットを参照)。これは 2 次元の曲線であるため、その次元は 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 の pp は多変量、つまり一変量 pp 関数のテンソル積とすることもできます。このような多変量 pp の pp 型は、breaks が各変数のブレーク シーケンスを含む cell 配列となり、coefs が多次元配列となるため、少しだけ複雑度が増します。このような非ランダム関数をゼロから構成するのはかなり困難であり、特にツールボックスの目的は pp 関数に関するある条件からその関数の作成を支援することであるため、ここでは行いません。たとえば、次の Figure の球体は 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')

pp 型のスプラインにより効率的に "評価" できますが、通常は、最初に "B 型"、つまり B スプラインの線形結合の表現を決定することで、あるデータからのスプラインの "作成" はより効率的に処理されます。詳細については、B 型の作成と操作を参照してください。