補間法（スプライン）

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double h[], b[], d[], g[], u[], r[], x[], y[], sp, x1, y1;
		int n;
					// 設定
		n = 5;
		h = new double [n];
		b = new double [n];
		d = new double [n];
		g = new double [n];
		u = new double [n];
		r = new double [n+1];
		x = new double [n+1];
		y = new double [n+1];

		x[0] = 0.0;
		x[1] = 1.0;
		x[2] = 2.0;
		x[3] = 3.0;
		x[4] = 4.0;
		x[5] = 5.0;
		y[0] = 0.0;
		y[1] = 1.0;
		y[2] = 4.0;
		y[3] = 9.0;
		y[4] = 16.0;
		y[5] = 25.0;
					// 実行と出力
		x1 = 0.0;
		sp = 0.1;
		while (x1 <= 5.01) {
			y1 = spline(n, x, y, x1, h, b, d, g, u, r);
			System.out.println(x1 + " " + y1);
			x1 += sp;
		}
	}

	/**********************************/
	/* ３次スプライン関数による補間   */
	/*      n : 区間の数              */
	/*      x, y : 点の座標           */
	/*      x1 : 補間値を求める値     */
	/*      h, b, d, g, u, r : 作業域 */
	/*      return : 補間値           */
	/*      coded by Y.Suganuma       */
	/**********************************/
	static double spline(int n, double x[], double y[], double x1,
           double h[], double b[], double d[], double g[], double u[], double r[])
	{
		int i = -1, i1, k;
		double y1, qi, si, xx;
					// 区間の決定
		for (i1 = 1; i1 < n && i < 0; i1++) {
			if (x1 < x[i1])
				i = i1 - 1;
		}
		if (i < 0)
			i = n - 1;
					// ステップ１
		for (i1 = 0; i1 < n; i1++)
			h[i1] = x[i1+1] - x[i1];
		for (i1 = 1; i1 < n; i1++) {
			b[i1] = 2.0 * (h[i1] + h[i1-1]);
			d[i1] = 3.0 * ((y[i1+1] - y[i1]) / h[i1] - (y[i1] - y[i1-1]) / h[i1-1]);
		}
					// ステップ２
		g[1] = h[1] / b[1];
		for (i1 = 2; i1 < n-1; i1++)
			g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1]);
		u[1] = d[1] / b[1];
		for (i1 = 2; i1 < n; i1++)
			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1]);
					// ステップ３
		k      = (i > 1) ? i : 1;
		r[0]   = 0.0;
		r[n]   = 0.0;
		r[n-1] = u[n-1];
		for (i1 = n-2; i1 >= k; i1--)
			r[i1] = u[i1] - g[i1] * r[i1+1];
					// ステップ４
		xx = x1 - x[i];
		qi = (y[i+1] - y[i]) / h[i] - h[i] * (r[i+1] + 2.0 * r[i]) / 3.0;
		si = (r[i+1] - r[i]) / (3.0 * h[i]);
		y1 = y[i] + xx * (qi + xx * (r[i]  + si * xx));

		return y1;
	}
}
		

以下は，C++のクラスを使用した例に対応したプログラムです．

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double x[], y[], sp, x1, y1;
		int n;
					// 設定
		n = 5;
		x = new double [n+1];
		y = new double [n+1];

		x[0] = 0.0;
		x[1] = 1.0;
		x[2] = 2.0;
		x[3] = 3.0;
		x[4] = 4.0;
		x[5] = 5.0;
		y[0] = 0.0;
		y[1] = 1.0;
		y[2] = 4.0;
		y[3] = 9.0;
		y[4] = 16.0;
		y[5] = 25.0;

		Spline spline = new Spline(n, x, y);
					// 実行と出力
		x1 = 0.0;
		sp = 0.1;
		while (x1 <= 5.01) {
			y1 = spline.value(x1);
			System.out.println(x1 + " " + y1);
			x1 += sp;
		}
	}
}

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/
class Spline
{
	private double h[], b[], d[], g[], u[], r[], q[], s[], x[], y[];
	private int n;

	/**************************/
	/* コンストラクタ         */
	/*      n1 : 区間の数     */
	/*      x1, y1 : 点の座標 */
	/**************************/
	Spline(int n1, double x1[], double y1[])
	{
		int i1;
						// 領域の確保
		n = n1;

		h = new double [n];
		b = new double [n];
		d = new double [n];
		g = new double [n];
		u = new double [n];
		q = new double [n];
		s = new double [n];
		r = new double [n+1];
		x = new double [n+1];
		y = new double [n+1];

		for (i1 = 0; i1 <= n; i1++) {
			x[i1] = x1[i1];
			y[i1] = y1[i1];
		}
						// ステップ１
		for (i1 = 0; i1 < n; i1++)
			h[i1] = x[i1+1] - x[i1];
		for (i1 = 1; i1 < n; i1++) {
			b[i1] = 2.0 * (h[i1] + h[i1-1]);
			d[i1] = 3.0 * ((y[i1+1] - y[i1]) / h[i1] - (y[i1] - y[i1-1]) / h[i1-1]);
		}
						// ステップ２
		g[1] = h[1] / b[1];
		for (i1 = 2; i1 < n-1; i1++)
			g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1]);
		u[1] = d[1] / b[1];
		for (i1 = 2; i1 < n; i1++)
			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1]);
						// ステップ３
		r[0]   = 0.0;
		r[n]   = 0.0;
		r[n-1] = u[n-1];
		for (i1 = n-2; i1 >= 1; i1--)
			r[i1] = u[i1] - g[i1] * r[i1+1];
						// ステップ４
		for (i1 = 0; i1 < n; i1++) {
			q[i1] = (y[i1+1] - y[i1]) / h[i1] - h[i1] * (r[i1+1] + 2.0 * r[i1]) / 3.0;
			s[i1] = (r[i1+1] - r[i1]) / (3.0 * h[i1]);
		}
	}

	/******************************/
	/* 補間値の計算               */
	/*      x1 : 補間値を求める値 */
	/*      return : 補間値       */
	/******************************/
	double value(double x1)
	{
		int i = -1, i1;
		double y1, xx;
						// 区間の決定
		for (i1 = 1; i1 < n && i < 0; i1++) {
			if (x1 < x[i1])
				i = i1 - 1;
		}
		if (i < 0)
			i = n - 1;
						// 計算
		xx = x1 - x[i];
		y1 = y[i] + xx * (q[i] + xx * (r[i]  + s[i] * xx));

		return y1;
	}
}