補間法(3次スプライン関数)

/********************************/
/* 補間法(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;
	}
}