静岡理工科大学 菅沼ホーム 目次 索引

補間法(スプライン)

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,3 次スプライン関数によってスプライン補間するためのものです.なお,このプログラムは,クラスを使用して記述しています( JavaScript の場合はオブジェクト).

  1. C++

    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      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;
    }
    
    /**********************************/
    /* 3次スプライン関数による補間   */
    /*      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;
    					// ステップ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]);
    	}
    					// ステップ2
    	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]);
    					// ステップ3
    	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];
    					// ステップ4
    	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;
    }
    			
    クラスを使用した場合
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    #include <stdio.h>
    
    /****************/
    /* クラスの定義 */
    /****************/
    class Spline
    {
    		double *r, *q, *s, *x, *y;
    		int n;
    
    	public:
    
    		/**************************/
    		/* コンストラクタ         */
    		/*      n1 : 区間の数     */
    		/*      x1, y1 : 点の座標 */
    		/**************************/
    		Spline(int n1, double *x1, double *y1)
    		{
    			int i1;
    						// 領域の確保
    			n = n1;
    
    			double *h = new double [n];
    			double *b = new double [n];
    			double *d = new double [n];
    			double *g = new double [n];
    			double *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];
    			}
    						// ステップ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]);
    			}
    						// ステップ2
    			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]);
    						// ステップ3
    			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];
    						// ステップ4
    			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]);
    			}
    
    			delete [] h;
    			delete [] b;
    			delete [] d;
    			delete [] g;
    			delete [] u;
    		}
    
    		/****************/
    		/* デストラクタ */
    		/****************/
    		~Spline()
    		{
    			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;
    }
    			

  2. Java

      アプレット版では,任意のデータに対して,スプライン補間法による計算結果を画面上で求めることができます.
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      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[], 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];
    		double x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    		double y[] = {0.0, 1.0, 4.0, 9.0, 16.0, 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;
    		}
    	}
    
    	/**********************************/
    	/* 3次スプライン関数による補間   */
    	/*      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;
    					// ステップ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]);
    		}
    					// ステップ2
    		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]);
    					// ステップ3
    		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];
    					// ステップ4
    		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;
    	}
    }
    			
    クラスを使用した場合
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    import java.io.*;
    
    public class Test_c {
    	public static void main(String args[]) throws IOException
    	{
    		double sp, x1, y1;
    		int n;
    					// 設定
    		n = 5;
    		double x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    		double y[] = {0.0, 1.0, 4.0, 9.0, 16.0, 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;
    		}
    	}
    }
    
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    class Spline
    {
    	private double r[], q[], s[], x[], y[];
    	private int n;
    
    	/**************************/
    	/* コンストラクタ         */
    	/*      n1 : 区間の数     */
    	/*      x1, y1 : 点の座標 */
    	/**************************/
    	Spline(int n1, double x1[], double y1[])
    	{
    		int i1;
    						// 領域の確保
    		n = n1;
    
    		double h[] = new double [n];
    		double b[] = new double [n];
    		double d[] = new double [n];
    		double g[] = new double [n];
    		double 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];
    		}
    						// ステップ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]);
    		}
    						// ステップ2
    		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]);
    						// ステップ3
    		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];
    						// ステップ4
    		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;
    	}
    }
    			
    アプレット版
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Spline extends Applet implements ActionListener {
    	TextField order, step;
    	TextArea ta1, ta2;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("区間の数(n):"));
    		order = new TextField("5", 3);
    		pn1.add(order);
    
    		pn1.add(new Label("ステップ:"));
    		step = new TextField("0.1", 3);
    		pn1.add(step);
    
    		pn1.add(new Label(" "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn1.add(bt);
    					// 中央のパネル(データ)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		pn2.add(new Label("データ(x y,(n+1)個):"));
    		ta1 = new TextArea("0.0 0.0\n1.0 1.0\n2.0 4.0\n3.0 9.0\n4.0 16.0\n5.0 25.0", 10, 30);
    		pn2.add(ta1);
    					// 下のパネル(補間値)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		ta2 = new TextArea(10, 30);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    					// 問題の設定
    			int n = Integer.parseInt(order.getText());
    			double sp = Double.parseDouble(step.getText());
    			double h[] = new double [n];
    			double b[] = new double [n];
    			double d[] = new double [n];
    			double g[] = new double [n];
    			double u[] = new double [n];
    			double r[] = new double [n+1];
    			double x[] = new double [n+1];
    			double y[] = new double [n+1];
    
    			String s = ta1.getText();
    			StringTokenizer str = new StringTokenizer(s, " \n");
    			int k = 0;
    			while (str.hasMoreTokens() && k < n+1) {
    				x[k] = Double.parseDouble(str.nextToken());
    				y[k] = Double.parseDouble(str.nextToken());
    				k++;
    			}
    					// 計算
    			ta2.setText("");
    			double x1 = x[0];
    			while (x1 < x[n]+0.01*sp) {
    				double y1 = spline(n, x, y, x1, h, b, d, g, u, r);
    				ta2.append("x " + String.format("%.5f",x1) + " y " + String.format("%.5f",y1) + "\n");
    				x1 += sp;
    			}
    		}
    	}
    
    	/**********************************/
    	/* 3次スプライン関数による補間   */
    	/*      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;
    					// ステップ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]);
    		}
    					// ステップ2
    		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]);
    					// ステップ3
    		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];
    					// ステップ4
    		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;
    	}
    }
    			

  3. JavaScript

      ここをクリックオブジェクトを使用した場合)すると,任意のデータに対して,スプライン補間法による計算結果を画面上で求めることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>補間法(3次スプライン関数)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		function main()
    		{
    					// データの設定
    			let n  = parseInt(document.getElementById("order").value);
    			let sp = parseFloat(document.getElementById("step").value);
    			let h  = new Array();
    			let b  = new Array();
    			let d  = new Array();
    			let g  = new Array();
    			let u  = new Array();
    			let r  = new Array();
    			let x  = new Array();
    			let y  = new Array();
    			let s = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < 2*(n+1); i1 += 2) {
    				x[k] = parseFloat(s[i1]);
    				y[k] = parseFloat(s[i1+1]);
    				k++;
    			}
    					// 出力
    			let str = "";
    			let x1 = x[0];
    			while (x1 < x[n]+0.01*sp) {
    				let y1 = spline(n, x, y, x1, h, b, d, g, u, r);
    				str = str + "x " + x1 + " y " + y1 + "\n";
    				x1 += sp;
    			}
    			document.getElementById("ans").value = str;
    		}
    
    		/**********************************/
    		/* 3次スプライン関数による補間   */
    		/*      n : 区間の数              */
    		/*      x, y : 点の座標           */
    		/*      x1 : 補間値を求める値     */
    		/*      h, b, d, g, u, r : 作業域 */
    		/*      return : 補間値           */
    		/*      coded by Y.Suganuma       */
    		/**********************************/
    		function spline(n, x, y, x1, h, b, d, g, u, r)
    		{
    			let i = -1;
    			let i1;
    			let k;
    			let y1;
    			let qi;
    			let si;
    			let xx;
    						// 区間の決定
    			for (i1 = 1; i1 < n && i < 0; i1++) {
    				if (x1 < x[i1])
    					i = i1 - 1;
    			}
    			if (i < 0)
    				i = n - 1;
    						// ステップ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]);
    			}
    						// ステップ2
    			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]);
    						// ステップ3
    			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];
    						// ステップ4
    			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;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>補間法(3次スプライン関数)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,6 点 ( 0, 0 ), ( 1, 1 ), ( 2, 4 ) ( 3, 9 ), ( 4, 16 ), ( 5, 25 ) を使用して,3 次スプライン関数による補間を実行するための値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		区間の数(n):<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="5"> 
    		ステップ:<INPUT ID="step" STYLE="font-size: 100%;" TYPE="text" SIZE="5" VALUE="0.1"> 
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		データ(x y,(n+1)個):<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">0.0 0.0
    1.0 1.0
    2.0 4.0
    3.0 9.0
    4.0 16.0
    5.0 25.0</TEXTAREA><BR><BR>
    		<TEXTAREA ID="ans" COLS="50" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			
    オブジェクトを使用した場合
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>補間法(3次スプライン関数)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		spline = null;
    
    		function main()
    		{
    					// データの設定
    			let n  = parseInt(document.getElementById("order").value);
    			let sp = parseFloat(document.getElementById("step").value);
    			let x  = new Array();
    			let y  = new Array();
    			let s = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < 2*(n+1); i1 += 2) {
    				x[k] = parseFloat(s[i1]);
    				y[k] = parseFloat(s[i1+1]);
    				k++;
    			}
    
    			spline = new Spline(n, x, y);
    					// 出力
    			let str = "";
    			let x1 = x[0];
    			while (x1 < x[n]+0.01*sp) {
    				let y1 = spline.value(x1);
    				str = str + "x " + x1 + " y " + y1 + "\n";
    				x1 += sp;
    			}
    			document.getElementById("ans").value = str;
    		}
    
    		/********************************/
    		/* 3次スプライン関数による補間 */
    		/*      n1 : 区間の数           */
    		/*      x1, y1 : 点の座標       */
    		/*      coded by Y.Suganuma     */
    		/********************************/
    		function Spline(n1, x1, y1)
    		{
    						// 領域の確保
    			let h = new Array();
    			let b = new Array();
    			let d = new Array();
    			let g = new Array();
    			let u = new Array();
    			this.n = n1;
    			this.q = new Array();
    			this.s = new Array();
    			this.r = new Array();
    			this.x = new Array();
    			this.y = new Array();
    			for (let i1 = 0; i1 <= this.n; i1++) {
    				this.x[i1] = x1[i1];
    				this.y[i1] = y1[i1];
    			}
    						// ステップ1
    			for (let i1 = 0; i1 < this.n; i1++)
    				h[i1] = this.x[i1+1] - this.x[i1];
    			for (let i1 = 1; i1 < this.n; i1++) {
    				b[i1] = 2.0 * (h[i1] + h[i1-1]);
    				d[i1] = 3.0 * ((this.y[i1+1] - this.y[i1]) / h[i1] - (this.y[i1] - this.y[i1-1]) / h[i1-1]);
    			}
    						// ステップ2
    			g[1] = h[1] / b[1];
    			for (let i1 = 2; i1 < this.n-1; i1++)
    				g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1]);
    			u[1] = d[1] / b[1];
    			for (let i1 = 2; i1 < this.n; i1++)
    				u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1]);
    						// ステップ3
    			this.r[0]        = 0.0;
    			this.r[this.n]   = 0.0;
    			this.r[this.n-1] = u[this.n-1];
    			for (let i1 = this.n-2; i1 >= 1; i1--)
    				this.r[i1] = u[i1] - g[i1] * this.r[i1+1];
    						// ステップ4
    			for (let i1 = 0; i1 < this.n; i1++) {
    				this.q [i1] = (this.y[i1+1] - this.y[i1]) / h[i1] - h[i1] * (this.r[i1+1] + 2.0 * this.r[i1]) / 3.0;
    				this.s[i1] = (this.r[i1+1] - this.r[i1]) / (3.0 * h[i1]);
    			}
    		}
    
    		/******************************/
    		/* 補間値の計算               */
    		/*      x1 : 補間値を求める値 */
    		/*      return : 補間値       */
    		/******************************/
    		Spline.prototype.value = function(x1)
    		{
    						// 区間の決定
    			let i = -1;
    			for (let i1 = 1; i1 < spline.n && i < 0; i1++) {
    				if (x1 < spline.x[i1])
    					i = i1 - 1;
    			}
    			if (i < 0)
    				i = spline.n - 1;
    						// 計算
    			let xx = x1 - spline.x[i];
    			let y1 = spline.y[i] + xx * (spline.q[i] + xx * (spline.r[i]  + spline.s[i] * xx));
    
    			return y1;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>補間法(3次スプライン関数)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,6 点 ( 0, 0 ), ( 1, 1 ), ( 2, 4 ) ( 3, 9 ), ( 4, 16 ), ( 5, 25 ) を使用して,3 次スプライン関数による補間を実行するための値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		区間の数(n):<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="5"> 
    		ステップ:<INPUT ID="step" STYLE="font-size: 100%;" TYPE="text" SIZE="5" VALUE="0.1"> 
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		データ(x y,(n+1)個):<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">0.0 0.0
    1.0 1.0
    2.0 4.0
    3.0 9.0
    4.0 16.0
    5.0 25.0</TEXTAREA><BR><BR>
    		<TEXTAREA ID="ans" COLS="50" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    
    					// 設定
    	$n = 5;
    	$h = array($n);
    	$b = array($n);
    	$d = array($n);
    	$g = array($n);
    	$u = array($n);
    	$r = array($n+1);
    	$x = array($n+1);
    	$y = array($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;
    	}
    
    /**********************************/
    /* 3次スプライン関数による補間   */
    /*      n : 区間の数              */
    /*      x, y : 点の座標           */
    /*      x1 : 補間値を求める値     */
    /*      h, b, d, g, u, r : 作業域 */
    /*      return : 補間値           */
    /*      coded by Y.Suganuma       */
    /**********************************/
    function spline($n, $x, $y, $x1, $h, $b, $d, $g, $u, $r)
    {
    	$i = -1;
    					// 区間の決定
    	for ($i1 = 1; $i1 < $n && $i < 0; $i1++) {
    		if ($x1 < $x[$i1])
    			$i = $i1 - 1;
    	}
    	if ($i < 0)
    		$i = $n - 1;
    					// ステップ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]);
    	}
    					// ステップ2
    	$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]);
    					// ステップ3
    	$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];
    					// ステップ4
    	$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;
    }
    
    ?>
    			
    クラスを使用した場合
    <?php
    
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    
    /****************/
    /* クラスの定義 */
    /****************/
    class Spline
    {
    	private $r = array();
    	private $q = array();
    	private $s = array();
    	private $x = array();
    	private $y = array();
    	private $n;
    
    	/**************************/
    	/* コンストラクタ         */
    	/*      n1 : 区間の数     */
    	/*      x1, y1 : 点の座標 */
    	/**************************/
    	function Spline($n1, $x1, $y1)
    	{
    					// 領域の確保
    		$this->n = $n1;
    
    		$h = array();
    		$b = array();
    		$d = array();
    		$g = array();
    		$u = array();
    
    		for ($i1 = 0; $i1 <= $this->n; $i1++) {
    			$this->x[$i1] = $x1[$i1];
    			$this->y[$i1] = $y1[$i1];
    		}
    					// ステップ1
    		for ($i1 = 0; $i1 < $this->n; $i1++)
    			$h[$i1] = $this->x[$i1+1] - $this->x[$i1];
    		for ($i1 = 1; $i1 < $this->n; $i1++) {
    			$b[$i1] = 2.0 * ($h[$i1] + $h[$i1-1]);
    			$d[$i1] = 3.0 * (($this->y[$i1+1] - $this->y[$i1]) / $h[$i1] - ($this->y[$i1] - $this->y[$i1-1]) / $h[$i1-1]);
    		}
    					// ステップ2
    		$g[1] = $h[1] / $b[1];
    		for ($i1 = 2; $i1 < $this->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 < $this->n; $i1++)
    			$u[$i1] = ($d[$i1] - $h[$i1-1] * $u[$i1-1]) / ($b[$i1] - $h[$i1-1] * $g[$i1-1]);
    					// ステップ3
    		$this->r[0]          = 0.0;
    		$this->r[$this->n]   = 0.0;
    		$this->r[$this->n-1] = $u[$this->n-1];
    		for ($i1 = $this->n-2; $i1 >= 1; $i1--)
    			$this->r[$i1] = $u[$i1] - $g[$i1] * $this->r[$i1+1];
    					// ステップ4
    		for ($i1 = 0; $i1 < $this->n; $i1++) {
    			$this->q[$i1] = ($this->y[$i1+1] - $this->y[$i1]) / $h[$i1] - $h[$i1] * ($this->r[$i1+1] + 2.0 * $this->r[$i1]) / 3.0;
    			$this->s[$i1] = ($this->r[$i1+1] - $this->r[$i1]) / (3.0 * $h[$i1]);
    		}
    	}
    
    	/******************************/
    	/* 補間値の計算               */
    	/*      x1 : 補間値を求める値 */
    	/*      return : 補間値       */
    	/******************************/
    	function value($x1)
    	{
    		$i = -1;
    					// 区間の決定
    		for ($i1 = 1; $i1 < $this->n && $i < 0; $i1++) {
    			if ($x1 < $this->x[$i1])
    				$i = $i1 - 1;
    		}
    		if ($i < 0)
    			$i = $this->n - 1;
    					// 計算
    		$xx = $x1 - $this->x[$i];
    		$y1 = $this->y[$i] + $xx * ($this->q[$i] + $xx * ($this->r[$i]  + $this->s[$i] * $xx));
    
    		return $y1;
    	}
    }
    
    					// 設定
    	$n = 5;
    	$x = array($n+1);
    	$y = array($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 = new 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;
    	}
    
    ?>
    			

  5. Ruby

    #*******************************/
    # 補間法(3次スプライン関数) */
    #      coded by Y.Suganuma     */
    #*******************************/
    
    #*********************************/
    # 3次スプライン関数による補間   */
    #      n : 区間の数              */
    #      x, y : 点の座標           */
    #      x1 : 補間値を求める値     */
    #      h, b, d, g, u, r : 作業域 */
    #      return : 補間値           */
    #      coded by Y.Suganuma       */
    #*********************************/
    def spline(n, x, y, x1, h, b, d, g, u, r)
    					# 区間の決定
    	i = -1
    	for i1 in 1 ... n
    		if x1 < x[i1]
    			i = i1 - 1
    		end
    		if i >= 0
    			break
    		end
    	end
    	if i < 0
    		i = n - 1
    	end
    					# ステップ1
    	for i1 in 0 ... n
    		h[i1] = x[i1+1] - x[i1]
    	end
    	for i1 in 1 ... n
    		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])
    	end
    					# ステップ2
    	g[1] = h[1] / b[1]
    	for i1 in 2 ... n-1
    		g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1])
    	end
    	u[1] = d[1] / b[1]
    	for i1 in 2 ... n
    		u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1])
    	end
    					# ステップ3
    	k      = (i > 1) ? i : 1
    	r[0]   = 0.0
    	r[n]   = 0.0
    	r[n-1] = u[n-1]
    	i1     = n - 2
    	while i1 >= k
    		r[i1] = u[i1] - g[i1] * r[i1+1]
    		i1 -= 1;
    	end
    					# ステップ4
    	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
    end
    
    					# 設定
    n = 5
    h = Array.new(n)
    b = Array.new(n)
    d = Array.new(n)
    g = Array.new(n)
    u = Array.new(n)
    r = Array.new(n+1)
    x = Array[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
    y = Array[0.0, 1.0, 4.0, 9.0, 16.0, 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
    end
    			
    クラスを使用した場合
    #*******************************/
    # 補間法(3次スプライン関数) */
    #      coded by Y.Suganuma     */
    #*******************************/
    
    #***************/
    # クラスの定義 */
    #***************/
    class Spline
    
    	#*************************/
    	# コンストラクタ         */
    	#      n1 : 区間の数     */
    	#      x1, y1 : 点の座標 */
    	#*************************/
    	def initialize(n1, x1, y1)
    					# 領域の確保
    		h = Array.new(n1)
    		b = Array.new(n1)
    		d = Array.new(n1)
    		g = Array.new(n1)
    		u = Array.new(n1)
    		@_n = n1
    		@_q = Array.new(@_n)
    		@_s = Array.new(@_n)
    		@_r = Array.new(@_n+1)
    		@_x = Array.new(@_n+1)
    		@_y = Array.new(@_n+1)
    
    		for i1 in 0 ... @_n+1
    			@_x[i1] = x1[i1]
    			@_y[i1] = y1[i1]
    		end
    					# ステップ1
    		for i1 in 0 ... @_n
    			h[i1] = @_x[i1+1] - @_x[i1]
    		end
    		for i1 in 1 ... @_n
    			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])
    		end
    					# ステップ2
    		g[1] = h[1] / b[1]
    		for i1 in 2 ... @_n-1
    			g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1])
    		end
    		u[1] = d[1] / b[1]
    		for i1 in 2 ... @_n
    			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1])
    		end
    					# ステップ3
    		@_r[0]     = 0.0
    		@_r[@_n]   = 0.0
    		@_r[@_n-1] = u[@_n-1]
    		i1         = @_n - 2
    		while i1 >= 1
    			@_r[i1] = u[i1] - g[i1] * @_r[i1+1]
    			i1     -= 1
    		end
    					# ステップ4
    		for i1 in 0 ... @_n
    			@_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])
    		end
    	end
    
    	#*****************************/
    	# 補間値の計算               */
    	#      x1 : 補間値を求める値 */
    	#      return : 補間値       */
    	#*****************************/
    	def value(x1)
    					# 区間の決定
    		i = -1
    		for i1 in 1 ... @_n
    			if (x1 < @_x[i1])
    				i = i1 - 1
    			end
    			if i >= 0
    				break
    			end
    		end
    		if (i < 0)
    			i = @_n - 1
    		end
    					# 計算
    		xx = x1 - @_x[i]
    		y1 = @_y[i] + xx * (@_q[i] + xx * (@_r[i]  + @_s[i] * xx))
    
    		return y1
    	end
    end
    
    					# 設定
    n = 5
    x = Array[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
    y = Array[0.0, 1.0, 4.0, 9.0, 16.0, 25.0]
    
    spline = Spline.new(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
    end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    #*********************************/
    # 3次スプライン関数による補間   */
    #      n : 区間の数              */
    #      x, y : 点の座標           */
    #      x1 : 補間値を求める値     */
    #      h, b, d, g, u, r : 作業域 */
    #      return : 補間値           */
    #      coded by Y.Suganuma       */
    #*********************************/
    def spline(n, x, y, x1, h, b, d, g, u, r) :
    					# 区間の決定
    	i = -1
    	for i1 in range(1, n) :
    		if x1 < x[i1] :
    			i = i1 - 1
    		if i >= 0 :
    			break
    	if i < 0 :
    		i = n - 1
    					# ステップ1
    	for i1 in range(0, n) :
    		h[i1] = x[i1+1] - x[i1]
    	for i1 in range(1, n) :
    		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])
    					# ステップ2
    	g[1] = h[1] / b[1]
    	for i1 in range(2, n-1) :
    		g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1])
    	u[1] = d[1] / b[1]
    	for i1 in range(2, n) :
    		u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1])
    					# ステップ3
    	if i1 > 1 :
    		k = i
    	else :
    		k = 1
    	r[0]   = 0.0
    	r[n]   = 0.0
    	r[n-1] = u[n-1]
    	for i1 in range(n-2, k-1, -1) :
    		r[i1] = u[i1] - g[i1] * r[i1+1]
    					# ステップ4
    	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
    
    #*******************************/
    # 補間法(3次スプライン関数) */
    #      coded by Y.Suganuma     */
    #*******************************/
    					# 設定
    n = 5
    x = np.array([0, 1, 2, 3, 4, 5], np.float)
    y = np.array([0, 1, 4, 9, 16, 25], np.float)
    h = np.empty(n, np.float)
    b = np.empty(n, np.float)
    d = np.empty(n, np.float)
    g = np.empty(n, np.float)
    u = np.empty(n, np.float)
    r = np.empty(n+1, np.float)
    					# 実行と出力
    x1 = 0.0
    sp = 0.1
    while x1 <= 5.01 :
    	y1 = spline(n, x, y, x1, h, b, d, g, u, r)
    	print(x1, y1)
    	x1 += sp
    			
    クラスを使用した場合
    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    ############################################
    # クラスSplineの定義
    #      coded by Y.Suganuma
    ############################################
    
    class Spline :
    
    	#########################
    	# コンストラクタ
    	#      n1 : 区間の数
    	#      x1, y1 : 点の座標
    	#########################
    
    	def __init__(self, n1, x1, y1) :
    			# 設定
    		h = np.empty(n1, np.float)
    		b = np.empty(n1, np.float)
    		d = np.empty(n1, np.float)
    		g = np.empty(n1, np.float)
    		u = np.empty(n1, np.float)
    		self.n = n1
    		self.q = np.empty(self.n, np.float)
    		self.s = np.empty(self.n, np.float)
    		self.r = np.empty(self.n+1, np.float)
    		self.x = np.empty(self.n+1, np.float)
    		self.y = np.empty(self.n+1, np.float)
    
    		for i1 in range(0, self.n+1) :
    			self.x[i1] = x1[i1]
    			self.y[i1] = y1[i1]
    			# ステップ1
    		for i1 in range(0, self.n) :
    			h[i1] = self.x[i1+1] - self.x[i1]
    		for i1 in range(1, self.n) :
    			b[i1] = 2.0 * (h[i1] + h[i1-1])
    			d[i1] = 3.0 * ((self.y[i1+1] - self.y[i1]) / h[i1] - (self.y[i1] - self.y[i1-1]) / h[i1-1])
    			# ステップ2
    		g[1] = h[1] / b[1]
    		for i1 in range(2, self.n-1) :
    			g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1])
    		u[1] = d[1] / b[1]
    		for i1 in range(2, self.n) :
    			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1])
    			# ステップ3
    		self.r[0]        = 0.0
    		self.r[self.n]   = 0.0
    		self.r[self.n-1] = u[self.n-1]
    		for i1 in range(self.n-2, 0, -1) :
    			self.r[i1] = u[i1] - g[i1] * self.r[i1+1]
    			# ステップ4
    		for i1 in range(0, self.n) :
    			self.q[i1] = (self.y[i1+1] - self.y[i1]) / h[i1] - h[i1] * (self.r[i1+1] + 2.0 * self.r[i1]) / 3.0
    			self.s[i1] = (self.r[i1+1] - self.r[i1]) / (3.0 * h[i1])
    
    	#########################
    	# 補間値の計算
    	#      x1 : 補間値を求める値
    	#      return : 補間値
    	#########################
    
    	def value(self, x1) :
    			# 区間の決定
    		i = -1
    		for i1 in range(1, self.n) :
    			if x1 < self.x[i1] :
    				i = i1 - 1
    				break
    		if i < 0 :
    			i = self.n - 1
    			# 計算
    		xx = x1 - self.x[i]
    		y1 = self.y[i] + xx * (self.q[i] + xx * (self.r[i]  + self.s[i] * xx))
    
    		return y1
    
    ############################################
    # 補間法(3次スプライン関数)
    #      coded by Y.Suganuma
    ############################################
    
    			# 設定
    n = 5
    x = np.array([0, 1, 2, 3, 4, 5], np.float)
    y = np.array([0, 1, 4, 9, 16, 25], np.float)
    
    spline = Spline(n, x, y)
    			# 実行と出力
    x1 = 0.0
    sp = 0.1
    while x1 <= 5.01 :
    	y1 = spline.value(x1)
    	print(x1, y1)
    	x1 += sp
    			

  7. C#

    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    					// 設定
    		int n = 5;
    		double[] h = new double [n];
    		double[] b = new double [n];
    		double[] d = new double [n];
    		double[] g = new double [n];
    		double[] u = new double [n];
    		double[] r = new double [n+1];
    		double[] x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    		double[] y = {0.0, 1.0, 4.0, 9.0, 16.0, 25.0};
    					// 実行と出力
    		double x1 = 0.0;
    		double sp = 0.1;
    		while (x1 <= 5.01) {
    			double y1 = spline(n, x, y, x1, h, b, d, g, u, r);
    			Console.WriteLine(x1 + " " + y1);
    			x1 += sp;
    		}
    	}
    
    	/**********************************/
    	/* 3次スプライン関数による補間   */
    	/*      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;
    		for (int i1 = 1; i1 < n && i < 0; i1++) {
    			if (x1 < x[i1])
    				i = i1 - 1;
    		}
    		if (i < 0)
    			i = n - 1;
    					// ステップ1
    		for (int i1 = 0; i1 < n; i1++)
    			h[i1] = x[i1+1] - x[i1];
    		for (int 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]);
    		}
    					// ステップ2
    		g[1] = h[1] / b[1];
    		for (int 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 (int i1 = 2; i1 < n; i1++)
    			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1]);
    					// ステップ3
    		int k  = (i > 1) ? i : 1;
    		r[0]   = 0.0;
    		r[n]   = 0.0;
    		r[n-1] = u[n-1];
    		for (int i1 = n-2; i1 >= k; i1--)
    			r[i1] = u[i1] - g[i1] * r[i1+1];
    					// ステップ4
    		double xx = x1 - x[i];
    		double qi = (y[i+1] - y[i]) / h[i] - h[i] * (r[i+1] + 2.0 * r[i]) / 3.0;
    		double si = (r[i+1] - r[i]) / (3.0 * h[i]);
    		double y1 = y[i] + xx * (qi + xx * (r[i]  + si * xx));
    
    		return y1;
    	}
    }
    			
    クラスを使用した場合
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    					// 設定
    		int n = 5;
    		double[] x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
    		double[] y = {0.0, 1.0, 4.0, 9.0, 16.0, 25.0};
    
    		Spline spline = new Spline(n, x, y);
    					// 実行と出力
    		double x1 = 0.0;
    		double sp = 0.1;
    		while (x1 <= 5.01) {
    			double y1 = spline.value(x1);
    			Console.WriteLine(x1 + " " + y1);
    			x1 += sp;
    		}
    	}
    }
    
    /********************************/
    /* 補間法(3次スプライン関数) */
    /*      coded by Y.Suganuma     */
    /********************************/
    class Spline
    {
    	private double[] q;
    	private double[] s;
    	private double[] r;
    	private double[] x;
    	private double[] y;
    	private int n;
    
    	/**************************/
    	/* コンストラクタ         */
    	/*      n1 : 区間の数     */
    	/*      x1, y1 : 点の座標 */
    	/**************************/
    	public Spline(int n1, double[] x1, double[] y1)
    	{
    						// 領域の確保
    		n = n1;
    
    		double[] h = new double [n];
    		double[] b = new double [n];
    		double[] d = new double [n];
    		double[] g = new double [n];
    		double[] 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 (int i1 = 0; i1 <= n; i1++) {
    			x[i1] = x1[i1];
    			y[i1] = y1[i1];
    		}
    						// ステップ1
    		for (int i1 = 0; i1 < n; i1++)
    			h[i1] = x[i1+1] - x[i1];
    		for (int 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]);
    		}
    						// ステップ2
    		g[1] = h[1] / b[1];
    		for (int 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 (int i1 = 2; i1 < n; i1++)
    			u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1]);
    						// ステップ3
    		r[0]   = 0.0;
    		r[n]   = 0.0;
    		r[n-1] = u[n-1];
    		for (int i1 = n-2; i1 >= 1; i1--)
    			r[i1] = u[i1] - g[i1] * r[i1+1];
    						// ステップ4
    		for (int 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 : 補間値       */
    	/******************************/
    	public double value(double x1)
    	{
    						// 区間の決定
    		int i = -1;
    		for (int i1 = 1; i1 < n && i < 0; i1++) {
    			if (x1 < x[i1])
    				i = i1 - 1;
    		}
    		if (i < 0)
    			i = n - 1;
    						// 計算
    		double xx = x1 - x[i];
    		double y1 = y[i] + xx * (q[i] + xx * (r[i]  + s[i] * xx));
    
    		return y1;
    	}
    }
    			

  8. VB

    '******************************'
    ' 補間法(3次スプライン関数) '
    '      coded by Y.Suganuma     '
    '******************************'
    Module Test
    
    	Sub Main()
    					' 設定
    		Dim n As Integer = 5
    		Dim h(n) As Double
    		Dim b(n) As Double
    		Dim d(n) As Double
    		Dim g(n) As Double
    		Dim u(n) As Double
    		Dim r(n+1) As Double
    		Dim x() As Double = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0}
    		Dim y() As Double = {0.0, 1.0, 4.0, 9.0, 16.0, 25.0}
    					' 実行と出力
    		Dim x1 As Double = 0.0
    		Dim sp As Double = 0.1
    		Do While x1 <= 5.01
    			Dim y1 As Double = spline(n, x, y, x1, h, b, d, g, u, r)
    			Console.WriteLine(x1 & " " & y1)
    			x1 += sp
    		Loop
    
    	End Sub
    
    	'********************************'
    	' 3次スプライン関数による補間   '
    	'      n : 区間の数              '
    	'      x, y : 点の座標           '
    	'      x1 : 補間値を求める値     '
    	'      h, b, d, g, u, r : 作業域 '
    	'      return : 補間値           '
    	'      coded by Y.Suganuma       '
    	'********************************'
    	Function spline(n As Integer, x() As Double, y() As Double, x1 As Double,
    	                h() As Double, b() As Double, d() As Double, g() As Double,
    	                u() As Double, r() As Double)
    
    					' 区間の決定
    		Dim i As Integer  = -1
    		Dim i1 As Integer = 1
    		Do While i1 < n and i < 0
    			If x1 < x(i1)
    				i = i1 - 1
    			End If
    			i1 += 1
    		Loop
    		If i < 0
    			i = n - 1
    		End If
    					' ステップ1
    		For i1 = 0 To n-1
    			h(i1) = x(i1+1) - x(i1)
    		Next
    		For i1 = 1 To n-1
    			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))
    		Next
    					' ステップ2
    		g(1) = h(1) / b(1)
    		For i1 = 2 To n-2
    			g(i1) = h(i1) / (b(i1) - h(i1-1) * g(i1-1))
    		Next
    		u(1) = d(1) / b(1)
    		For i1 = 2 To n-1
    			u(i1) = (d(i1) - h(i1-1) * u(i1-1)) / (b(i1) - h(i1-1) * g(i1-1))
    		Next
    					' ステップ3
    		Dim k As Integer = i
    		If i <= 1
    			k = 1
    		End If
    		r(0)   = 0.0
    		r(n)   = 0.0
    		r(n-1) = u(n-1)
    		For i1 = n-2 To k Step -1
    			r(i1) = u(i1) - g(i1) * r(i1+1)
    		Next
    					' ステップ4
    		Dim xx As Double = x1 - x(i)
    		Dim qi As Double = (y(i+1) - y(i)) / h(i) - h(i) * (r(i+1) + 2.0 * r(i)) / 3.0
    		Dim si As Double = (r(i+1) - r(i)) / (3.0 * h(i))
    		Dim y1 As Double = y(i) + xx * (qi + xx * (r(i)  + si * xx))
    
    		Return y1
    
    	End Function
    
    End Module
    			
    クラスを使用した場合
    '******************************'
    ' 補間法(3次スプライン関数) '
    '      coded by Y.Suganuma     '
    '******************************'
    Module Test
    
    	Sub Main()
    					' 設定
    		Dim n As Integer  = 5
    		Dim x() As Double = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0}
    		Dim y() As Double = {0.0, 1.0, 4.0, 9.0, 16.0, 25.0}
    	
    		Dim sp As Spline = new Spline(n, x, y)
    					' 実行と出力
    		Dim x1 As Double  = 0.0
    		Dim sp1 As Double = 0.1
    		Do While x1 <= 5.01
    			Dim y1 As Double = sp.value(x1)
    			Console.WriteLine(x1 & " " & y1)
    			x1 += sp1
    		Loop
    	End Sub
    
    	'******************************'
    	' 補間法(3次スプライン関数) '
    	'      coded by Y.Suganuma     '
    	'******************************'
    	Class Spline
    
    		Private n As Integer
    		Private q() As Double
    		Private s() As Double
    		Private r() As Double
    		Private x() As Double
    		Private y() As Double
    	
    		'************************'
    		' コンストラクタ         '
    		'      n1 : 区間の数     '
    		'      x1, y1 : 点の座標 '
    		'************************'
    		Public Sub New (n1 As Integer, x1() As Double, y1() As Double)
    							' 領域の確保
    			n = n1
    			ReDim q(n)
    			ReDim s(n)
    			ReDim r(n+1)
    			ReDim x(n+1)
    			ReDim y(n+1)
    	
    			Dim h(n) As Double
    			Dim b(n) As Double
    			Dim d(n) As Double
    			Dim g(n) As Double
    			Dim u(n) As Double
    	
    			For i1 As Integer = 0 To n
    				x(i1) = x1(i1)
    				y(i1) = y1(i1)
    			Next
     							' ステップ1
    			For i1 As Integer = 0 To n-1
    				h(i1) = x(i1+1) - x(i1)
    			Next
    			For i1 As Integer = 1 To n-1
    				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))
    			Next
    							' ステップ2
    			g(1) = h(1) / b(1)
    			For i1 As Integer = 2 To n-2
    				g(i1) = h(i1) / (b(i1) - h(i1-1) * g(i1-1))
    			Next
    			u(1) = d(1) / b(1)
    			For i1 As Integer = 2 To n-1
    				u(i1) = (d(i1) - h(i1-1) * u(i1-1)) / (b(i1) - h(i1-1) * g(i1-1))
    			Next
    							' ステップ3
    			r(0)   = 0.0
    			r(n)   = 0.0
    			r(n-1) = u(n-1)
    			For i1 As Integer = n-2 To 1 Step -1
    				r(i1) = u(i1) - g(i1) * r(i1+1)
    			Next
    							' ステップ4
    			For i1 As Integer = 0 To n-1
    				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))
    			Next
    
    		End Sub
    	
    		'****************************'
    		' 補間値の計算               '
    		'      x1 : 補間値を求める値 '
    		'      return : 補間値       '
    		'****************************'
    		Function value(x1 As Double) As Double
    							' 区間の決定
    			Dim i As Integer  = -1
    			Dim i1 As Integer = 1
    			Do While i1 < n and i < 0
    				If x1 < x(i1)
    					i = i1 - 1
    				End If
    				I1 += 1
    			Loop
    			If i < 0
    				i = n - 1
    			End If
    							' 計算
    			Dim xx As Double = x1 - x(i)
    			Dim y1 As Double = y(i) + xx * (q(i) + xx * (r(i)  + s(i) * xx))
    	
    			Return y1
    	
    		End Function
    	
    	End Class
    
    End Module
    			

静岡理工科大学 菅沼ホーム 目次 索引