deBoorの公式

上の記事で紹介したMATLABのcsapiはdeBoorの公式を用いて双3次スプラインを行っています。
この記事ではそのdeBoorの公式を自作してみたいと思います。

アルゴリズム

さて、ここでどのようにfxやfyを計算するかというと、
この記事の方法を用います。スプライン関数はx,y 変数ともに3次関数なので、この3次近似の方法を用いれば十分です。

M = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1;
 1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0;
 1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0;
 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
 0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3;
 0,0,0,0,1,0,0,0,2,0,0,0,3,0,0,0;
 0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0;
 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
 0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3;
 0,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0;
 0,1,2,3,0,0,0,0,0,0,0,0,0,0,0,0;
 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
 0,0,0,0,0,1,2,3,0,2,4,6,0,3,6,9;
 0,0,0,0,0,1,0,0,0,0,2,0,0,0,3,0;
 0,0,0,0,0,1,2,3,0,0,0,0,0,0,0,0;
 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0];

b = [f_Z(i,j);
     f_Z(i,j+1);
     f_Z(i+1,j);
     f_Z(i+1,j+1);
     -h*fx(i,j);
     -h*fx(i,j+1);
     -h*fx(i+1,j);
     -h*fx(i+1,j+1);
     -k*fy(i,j);
     -k*fy(i,j+1);
     -k*fy(i+1,j);
     -k*fy(i+1,j+1);
     h*k*fxy(i,j);
     h*k*fxy(i,j+1);
     h*k*fxy(i+1,j);
     h*k*fxy(i+1,j+1)];

(投げやりですが)なる行列M、ベクトルbについて、連立一次方程式Ma = bをとけばaが求まります。ここでh = x_i+1 - x_i, k = y_i+1 - y_i, a = (a_00, a_01, a_02, a_03, a_10, a_11, ...)となります。


前回の記事のpeaks関数をdeBoorの公式で双3次スプライン補間して、等高線グラフで比較したものです。ソースコードは下にあります。

ということなんですが、
matlabのcsapiを比較すると、少しだけずれているんですよね...
端点条件が違うのかなー

ソースコード