補間法（スプライン）

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/
#include <stdio.h>

double spline(int, double *, double *, double,
              double *, double *, double *, double *, double *, double *);

int main()
{
	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);
		printf("%f %f\n", x1, y1);
		x1 += sp;
	}

	delete [] h;
	delete [] b;
	delete [] d;
	delete [] g;
	delete [] u;
	delete [] r;
	delete [] x;
	delete [] y;

	return 0;
}

/**********************************/
/* ３次スプライン関数による補間   */
/*      n : 区間の数              */
/*      x, y : 点の座標           */
/*      x1 : 補間値を求める値     */
/*      h, b, d, g, u, r : 作業域 */
/*      return : 補間値           */
/*      coded by Y.Suganuma       */
/**********************************/
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;
}
		

以下は，クラスを使用して書いた例です．

/********************************/
/* 補間法（３次スプライン関数） */
/*      coded by Y.Suganuma     */
/********************************/
#include <stdio.h>

/****************/
/* クラスの定義 */
/****************/
class Spline
{
		double *h, *b, *d, *g, *u, *r, *q, *s, *x, *y;
		int n;

	public:

		/**************************/
		/* コンストラクタ         */
		/*      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]);
			}
		}

		/****************/
		/* デストラクタ */
		/****************/
		~Spline()
		{
			delete [] h;
			delete [] b;
			delete [] d;
			delete [] g;
			delete [] u;
			delete [] q;
			delete [] s;
			delete [] r;
			delete [] x;
			delete [] y;
		}

		/******************************/
		/* 補間値の計算               */
		/*      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;
		}
};

/*************/
/* main 関数 */
/*************/
int main()
{
	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(n, x, y);
					// 実行と出力
	x1 = 0.0;
	sp = 0.1;
	while (x1 <= 5.01) {
		y1 = spline.value(x1);
		printf("%f %f\n", x1, y1);
		x1 += sp;
	}

	delete [] x;
	delete [] y;

	return 0;
}