Main Content

3 次スプライン内挿

この例では、Curve Fitting Toolbox™ の csapi コマンドおよび csape コマンドを使用して、3 次スプライン内挿を作成する方法を示します。

CSAPI コマンド

以下のコマンドを入力します。

values = csapi(x,y,xx)

は、節点なしの端点条件を使用して、与えられたデータ (x,y) に対する 3 次スプライン内挿の xx における値を返します。この内挿は、ブレーク シーケンス x をもつ区分的三次関数であり、その 3 次部分が結合して 2 つの連続状態変数の微係数をもつ関数を形成します。"節点なし" の端点条件とは、最初と最後の内部ブレークで、3 階微分でも連続である (丸め誤差を除く) ことを意味します。

2 つのデータ点のみを指定すると、直線内挿になります。

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Figure contains an axes object. The axes object with title Interpolant to Two Points contains 2 objects of type line. One or more of the lines displays its values using only markers

3 つのデータ点を指定すると、放物線が得られます。

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

Figure contains an axes object. The axes object with title Interpolant to Three Points contains 2 objects of type line. One or more of the lines displays its values using only markers

より一般的な例として、4 つ以上のデータ点を指定すると、3 次スプラインが得られます。

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

Figure contains an axes object. The axes object with title Cubic Spline Interpolant to Five Points contains 2 objects of type line. One or more of the lines displays its values using only markers

CSAPI の出力を確認する方法

これらの内挿は適切であるように見えますが、csapi の処理が説明どおりであるかをどのように確認すれば良いでしょうか。

データ点をプロットしたところ、内挿がこれらの点を正しく通っていたことから、csapi の内挿は既に確認できています。ただし、3 次スプラインが得られたことを確認するには、想定される種類の 3 次スプラインのデータから始め、csapi でその 3 次スプラインが "再現" されるかどうか、つまり、データの取得元である 3 次スプラインに戻るかどうかを確認するのが最善です。

例: 切断べき関数

3 次スプライン関数を確認する簡単な例としては、打ち切られた 3 のべき乗、つまり、次の関数があります。

f(x)=((x-xi)+)3,

ここで、xi はブレークの 1 つであり、添字の "+" はコマンド subplus で与えられる "切断関数" を示します。

help subplus
 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

特定の選択肢 xi = 2 について、打ち切られた 3 のべき乗を次にプロットします。予想どおり、2 の左側は 0 であり、2 の右側は (x-2)^3 のように上昇しています。

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Figure contains an axes object. The axes object contains an object of type line.

ここで、この特定の 3 次スプラインをデータ サイト 0:6 で内挿し、スプラインの上に内挿を黒でプロットします。

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Interpolant to ((x-2)SubScript +) SuperScript 3 baseline contains 3 objects of type line. One or more of the lines displays its values using only markers

内挿誤差

2 つの関数を比較する場合、通常はそれらの差異をプロットする方がはるかに多くの情報が得られます。

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Error in Cubic Spline Interpolation to ((x-2)SubScript +) SuperScript 3 baseline contains an object of type line.

この差異の大きさを考慮に入れるため、最大データ値も計算します。これにより、この誤差がやむを得ない丸め誤差よりも悪くないことが示されます。

max_y = max(abs(y))
max_y = 64

打ち切られたべき乗を再現できない場合

追加のテストとして、csapi により生成されたサイト 0:6 の内挿が一致しない場合の打ち切られたべき乗を内挿します。たとえば、csapi は "節点なし" 条件を使用するため、内挿スプラインの最初の内部ブレークは、実は節点ではありません。そのため、この内挿はこのサイトで連続状態変数の微係数を 3 つもつことになります。つまり、その 3 階微分がそのサイトの周りでは不連続であるため、そのサイトを中心とする、打ち切られた 3 のべき乗を再現できないということになります。

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Figure contains an axes object. The axes object with title Error in Not-a-Knot blank Interpolant blank to blank ((x-1)SubScript +) SuperScript 3 baseline contains an object of type line.

1 は最初の内部節点であるため、この内挿ではアクティブではなくなります。

差の大きさは .18 ほどですが、1 から遠ざかるに従って減衰しています。これは、"3 次スプライン内挿が基本的に局所的" であることを説明しています。

値ではなく pp 型の使用

後での評価に適した形式、微分の計算に適した形式またはその他の操作に適した形式で、内挿 3 次スプラインを保持することができます。これは、次のような形式で csapi を呼び出すことで行います。

pp = csapi(x,y)

これは、内挿の pp 型を返します。この形式は、新しい点 xx で、次のコマンドによって評価できます。

values = fnval(pp,xx)

この内挿は、次のコマンドによって微分できます。

dpp = fnder(pp)

あるいは、次のコマンドで積分できます。

ipp = fnint(pp)

これらはそれぞれ、微分または積分の pp 型を返します。

例: 内挿の微分と積分

内挿の微分を示すため、次の打ち切られたべき乗の微分

f2(x)=3((x-2)+)2,

をプロットします (ここでも黄色)。次に、その上に、元の 3 乗の切断べき関数に対する内挿の微分をプロットします (ここでも黒)。

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Derivative of Interpolant to ((x-2)SubScript +) SuperScript 3 baseline contains 2 objects of type line.

繰り返しになりますが、差異をプロットする方が比較で多くの情報が得られます。また前回同様、これは丸め誤差と同じくらいの大きさです。

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Error in Derivative of interpolant to ((x-2)SubScript +) SuperScript 3 baseline contains an object of type line.

打ち切られたべき乗の 2 階微分は次のとおりです。

f2(x)=6(x-2)+

この関数と、元の関数に対する内挿の 2 階微分との間の差のプロットから、今度はジャンプがありますが、丸め誤差内に収まっていることが示されています。

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Error in Second Derivative of Interpolant to ((x-2)SubScript +) SuperScript 3 baseline contains an object of type line.

打ち切られたべき乗の積分は次のとおりです。

F2(x)=((x-2)+)4/4.

この関数と、元の関数に対する内挿の積分との間の差のプロットから、やはり誤差が丸め誤差内に収まっていることが示されています。

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Error in Integral of Interpolant to ((x-2)SubScript +) SuperScript 3 baseline contains an object of type line.

CSAPE コマンド

csapi と同様、csape コマンドは任意のデータに対する 3 次スプライン内挿を提供します。ただし、追加でさまざまな端点条件を指定できます。その最も単純なバージョンである

pp = csape(x,y)

はラグランジュ端点条件を使用します。これは、csapi で用いられる節点なし条件の代替手段によく使用されます。csape は内挿の値を直接返さず、pp 型のみを返します。

例として、次の関数に対する内挿を再度検討します。

f1(x)=((x-1)+)3,

これは、csapi でうまく再現されません。csape から求められる内挿の誤差 (赤) と共に、csapi によって返される節点なし内挿の誤差 (黒) をプロットします。

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

Figure contains an axes object. The axes object with title Error in Not-a-Knot vs. Lagrange End Conditions contains 2 objects of type line. These objects represent Not-a-Knot, Lagrange.

この場合、2 つの内挿にあまり差はありません。

その他の端点条件: '自然な' スプライン内挿

また、csape コマンドは、内挿 3 次スプラインのその他の種類の端点条件をいくつか指定する方法も提供します。たとえば、次のコマンド

pp = csape(x,y,'variational')

は、いわゆる '自然な' 端点条件を使用します。つまり、両端のブレークポイントで 2 階微分が 0 になります。

この手順では、'自然な' 3 次スプライン内挿を次の関数に適用し、

f2(x)=((x-2)+)3,

誤差をプロットする方法を示します。次のコードは、'variational' 引数に対応する別の引数構文を使用して '自然な' スプライン内挿を計算します。'second' により、csape が極値データ サイトで 2 階微分を既定値の 0 に設定することを指定します。

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title Error in 'Natural' blank Spline blank Interpolation blank to blank ((x-2)SubScript +) SuperScript 3 baseline contains an object of type line.

右端付近にある大きな誤差に注意してください。これは、'自然な' 端点条件がその位置でゼロの 2 階微分をもつことを暗黙的に示していることに起因します。

その他の端点条件: 2 階微分の指定

適切な 2 階微分を明示的に使用すると、誤差を小さくすることもできます。まず、打ち切られたべき乗の端点における適切な 2 階微分を計算します。

endcond = 6*subplus(x([1 end])-2);

次に、端点の 2 階微分が、計算した 2 階微分値と一致することを指定して、内挿を作成します。これは、データ値と共に、左の端点条件に endcond(1) を指定し、右に endcond(2) を指定して行います。

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Figure contains an axes object. The axes object with title Error in Spline Interpolation to ((x-1)SubScript +) SuperScript 3 baseline When Matching the 2nd Derivative at Ends contains an object of type line.

その他の端点条件: 勾配の指定

また、csape では端点の "勾配" を指定することもできます。これを "固定" (または "完全") 3 次スプライン内挿といいます。ステートメント

pp = csape(x,[sl,y,sr],'clamped')

によって、さらに左端のデータ サイトに勾配 sl、右端のデータ サイトに勾配 sr をもつデータ (x, y) への 3 次スプライン内挿を作成します。

その他の端点条件: 端点条件の混在

これまでの条件を混在させることも可能です。たとえば、前述の切断べき関数

f1(x)=((x-1)+)3

では、x=0 で勾配が 0、x=6 で 2 階微分が 30 となります (最後のデータ サイト)。

そのため、左端の勾配と右端の曲率を一致させることで、結果の内挿に誤差がなくなることが予想されます。

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Figure contains an axes object. The axes object with title Error in Spline Interpolation to ((x-1)SubScript +) SuperScript 3 baseline blank blank blank blank blank blank blank blank blank with blank Mixed blank End blank Conditions. contains an object of type line.

その他の端点条件: 周期的条件

"周期的な" 端点条件を指定することも可能です。たとえば、正弦関数の周期は 2*pi であり、サイト (pi/2)*(-2:2) における値は [0 -1 0 1 0] となります。正弦関数と、これらのサイトにおける周期的 3 次スプライン内挿との差異は 2 パーセントに過ぎません。これはまずまずの値です。

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

Figure contains an axes object. The axes object with title Error in Periodic Cubic Spline Interpolation to sin(x) contains a line object which displays its values using only markers.

CSAPI または CSAPE で明示的にカバーされない端点条件

csape の既定周辺条件を用いた内挿を作成してから、これにゼロ値および一部の周辺条件への内挿の適切なスカラー倍を加算することで、csapi または csape で明示的にカバーされない端点条件を処理できます。2 つの '非標準' 周辺条件を満たす必要がある場合は、まず 2 行 2 列の線形システムを解かなければならない可能性があります。

たとえば、次のデータに対する 3 次スプライン内挿 s を計算するとします。

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

次の条件を適用します。

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

これは点 e において s の 1 階微分および 2 階微分に対して行います。

このデータは、次の固有パラメーターをもつ周辺条件を満たす 4 次多項式から生成されます。

e = x(1);
a = 2; b = -3; c = 4;

この固有条件を満たす内挿を作成するには、最初に作成する内挿が既定端点条件

pp1 = csape(x,y);

およびその最初の多項式区分の 1 階微分をもつようにします。

dp1 = fnder(fnbrk(pp1,1));

さらに、e において勾配が 1 となることを指定し、ゼロ データ値に対する 3 次スプライン内挿を作成します。

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

また、その最初の多項式区分の 1 階微分も作成します。

dp0 = fnder(fnbrk(pp0,1));

次に、pp1 および pp0 の両方について lambda を計算します。

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

pp1 および pp0 の適切な線形結合を作成し、次の 3 次スプラインを求めます。

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

これは目的の条件のほか、右の端点における既定の端点条件も満たしています。この線形結合は、fncmb を使用して形成します。

s = fncmb(pp0,(c-lam1)/lam0,pp1);

内挿誤差のプロットは、e 付近では、s による 4 次多項式の近似の方が、既定の条件を使用した内挿 pp1 よりも優れていることを示しています。

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Default conditions, Nonstandard conditions.

次の条件

mu(s) := (D^3 s)(3) = 14.6

を内挿の 3 階微分 (4 次式はこの条件を満たす) に対して適用する場合、左の端点における 1 階微分がゼロであり、そのため pp0 から確実に独立している状態で、ゼロ値に対する 3 次スプライン内挿を追加作成します。

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

次に、次の線形結合の係数 d0 および d2 を求め、

s := pp1 + d0*pp0 + d2*pp2

線形システムを解きます。

lambda(s) = c

mu(s) = 14.6

pp0 および pp2 はどちらも、すべての内挿サイトで 0 になるため、s は選択したどの d0 および d2 についても任意のデータと一致することに注意してください。

試しに、MATLAB® エンコード機能で j=0:2 の場合の lambda(pp_j) および mu(pp_j) を計算するループを記述します。

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

pp0pp1 および pp2lambdamu の値が得られたら、適切な線形結合を定義する係数の値を求めます。

d = dd(:,[1,3])\([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

Figure contains an axes object. The axes object with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains a line object which displays its values using only markers.

再確認のため、この関数に対する完全 3 次スプライン内挿で得られたものを使用してこの誤差を計算します。

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

Figure contains an axes object. The axes object with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Nonstandard conditions, Endslope conditions.

この誤差は、終点付近でのみ (やや) 相違があり、pp0 および pp2 の両方がそれぞれの終点付近でのみ大きいという事実を示しています。

最終チェックとして、s が 3 における必要な 3 階微分条件を満たしていることを確認します。

fnval(fnder(s,3),3)
ans = 14.6000