3次スプライン補間
// 3次スプライン補間 | |
// 参考:高橋大輔,数値計算,岩波書店,1996. | |
// 入力 N+1個の点列 points = [[x0,y0],[x1,x2],...,[xn,yn]] | |
// 出力 各データ点間の区間における多項式の係数を含む2次元配列 | |
// x:x座標値で昇順ソートされた点列のx座標 | |
// y:x座標値で昇順ソートされた点列のy座標 | |
// a:区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1 | |
// b:区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1 | |
// c:区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1 | |
// d:区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1 | |
// ※三次式 y = a*x^3 + b*x^2 + c*x +d | |
function spline3Interpolation(points) { | |
// 入力点数が3未満ではスプライン曲線を定義できないので終了 | |
if(points.length<3) { | |
return null; | |
} | |
// x座標の値で昇順にソート | |
var tmp=numeric.clone(points); | |
var sortedPoints=tmp.sort( | |
function (a, b) { | |
if(a[0]<b[0]) { | |
return -1; | |
} else { | |
return 1; | |
} | |
} | |
); | |
// データ点 | |
var x=new Array(sortedPoints.length); | |
var y=new Array(sortedPoints.length); | |
// 曲線が定義される区間の数 | |
var N=points.length-1; | |
// 連立一次方程式を構成する行列 | |
var h=new Array(N); | |
var H=numeric.rep([N-1, N-1], 0); | |
var u=new Array(N+1); | |
var v=new Array(N-1); | |
// 各区間の係数 y = a*x^3 + b*x^2 + c*x +d | |
var a=new Array(N); | |
var b=new Array(N); | |
var c=new Array(N); | |
var d=new Array(N); | |
for(var i=0; i<sortedPoints.length; ++i) { | |
x[i]=sortedPoints[i][0]; | |
y[i]=sortedPoints[i][1]; | |
} | |
for(var i=0; i<N; ++i) { | |
h[i] = x[i+1] - x[i]; | |
} | |
H[0][0]=2*(h[0]+h[1]); | |
H[0][1]=h[1]; | |
for(var i=1; i<N-2; ++i) { | |
H[i][i-1]=h[i]; | |
H[i][i]=2*(h[i]+h[i+1]); | |
H[i][i+1]=h[i+1]; | |
} | |
H[N-2][N-3]=h[N-2]; | |
H[N-2][N-2]=2*(h[N-2]+h[N-1]); | |
for(var i=1; i<N; ++i) { | |
v[i-1]=6*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]); | |
} | |
var uBlock=numeric.solve(H, v); // from u[1] to u[N-1] | |
u[0]=0; | |
u[N]=0; | |
for(var i=0; i<N-1; ++i) { | |
u[i+1] = uBlock[i]; | |
} | |
for(var i=0; i<N; ++i) { | |
a[i]=(u[i+1]-u[i])/(6*(x[i+1]-x[i])); | |
} | |
for(var i=0; i<N; ++i) { | |
b[i]=0.5*u[i]; | |
} | |
for(var i=0; i<N; ++i) { | |
c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])-(x[i+1]-x[i])*(2*u[i]+u[i+1])/6; | |
} | |
for(var i=0; i<N; ++i) { | |
d[i]=y[i]; | |
} | |
return [x,y,a,b,c,d]; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment