上の記事で紹介した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を比較すると、少しだけずれているんですよね...
端点条件が違うのかなー