Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.