interpolation.java

数值分析算法源码(java) 这个学期一边学习java一边学习数值分析,因此用java写了一个数值分析算法的软件包numericalAnalysis. [说明] 适合使用者:会java的

package numericalAnalysis.interpolation;

import java.text.DecimalFormat;

/**
 * 本类定义了一个插值问题,包括多种插值公式和方法.对照书:数值分析第二版,史万明等编,北京理工大学出版社
 * 
 * @author 山
 * 
 */
public class Interpolation {
	private int n;
	private double[] x, y;
	private DecimalFormat fmt = new DecimalFormat("0.######");

	/**
	 * 新建一个插值问题.在外部类中定义n,赋值n=节点数减1,定义数组x[n+1]用于保存x[i]=xi(i=0,1,...,n);数组y同理
	 * 
	 * @param n
	 *            本参数是节点数减1的值
	 * @param x
	 *            x[i]=xi(i=0,1,...,n)
	 * @param y
	 *            y[i]=yi(i=0,1,...,n)
	 */
	public Interpolation(int n, double[] x, double[] y) {
		this.n = n;
		this.x = new double[this.n + 1];
		this.y = new double[this.n + 1];
		for (int i = 0; i < this.n + 1; i++)
			this.x[i] = x[i];
		for (int i = 0; i < this.n + 1; i++)
			this.y[i] = y[i];
	}

	/**
	 * 输出节点用以确认
	 * 
	 * @return 节点的xi,yi值
	 */
	public String print() {
		String output = "\n";

		output += "\n节点的xi:\n";
		for (int i = 0; i < n + 1; i++)
			output += "x" + i + "=(" + x[i] + ")   ";
		output += "\n节点的yi=f(xi):\n";
		for (int i = 0; i < n + 1; i++)
			output += "y" + i + "=f(x" + i + ")=(" + y[i] + ")   ";

		return output;
	}

	/**
	 * 插值求解f(x)的近似值Pn(x).用牛顿基本茶商公式的方法
	 * 
	 * @param x
	 *            y=Pn(x)的x
	 * @return y=Pn(x)的y
	 */
	public double getSolution(double x) {
		double y = 0;

		// compute y
		for (int i = 0; i < n + 1; i++) {
			double product = 1;
			for (int j = 0; j < i; j++)
				product *= x - this.x[j];
			y += product * differenceQuotient(i, this.x, this.y);
		}

		return y;
	}

	/**
	 * 计算一个区间的差商.定义n,赋值n=节点数减1,定义数组x[n+1]用于保存x[i]=xi(i=0,1,...,n);数组y同理
	 * 
	 * @param n
	 *            节点数减1的值
	 * @param x
	 *            x[i]=xi(i=0,1,...,n)
	 * @param y
	 *            y[i]=yi(i=0,1,...,n)
	 * @return [x0,xn]区间的差商
	 */
	public static double differenceQuotient(int n, double[] x, double[] y) {
		double result;
		if (n == 0)
			result = y[0];
		else {
			double[] x1 = new double[n], y1 = new double[n], x2 = new double[n], y2 = new double[n];

			for (int i = 0; i < n; i++) {
				x1[i] = x[i];
				y1[i] = y[i];
			}
			for (int i = 1; i < n + 1; i++) {
				x2[i - 1] = x[i];
				y2[i - 1] = y[i];
			}

			result = (differenceQuotient(n - 1, x2, y2) - differenceQuotient(
					n - 1, x1, y1))
					/ (x[n] - x[0]);
		}

		return result;
	}

	/**
	 * 输出差商表
	 * 
	 * @return 差商表
	 */
	public String differenceQuotientChart() {
		String output = "\n";

		output += "x\tf(x)\t";
		for (int i = 1; i < n + 1; i++) {
			output += "f[xi,";
			for (int j = 1; j < i + 1; j++) {
				output += "x(i+" + j + ")";
				if (j != i)
					output += ",";
			}
			output += "]\t";
		}
		output += "\n";

		for (int i = 0; i < n + 1; i++) {
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\t";
			for (int j = 1; j < i + 1; j++) {
				double[] tempX = new double[j + 1], tempY = new double[j + 1];
				for (int k = i - j, m = 0; k < i + 1; k++, m++) {
					tempX[m] = x[k];
					tempY[m] = y[k];
				}
				output += fmt.format(differenceQuotient(j, tempX, tempY))
						+ "\t";
			}
			output += "\n";
		}

		return output;
	}

	/**
	 * 牛顿基本差商公式
	 * 
	 * @param x
	 *            Pn(x)的x
	 * @return 结果
	 */
	public String newtonBasicDifferenceQuotient(double x) {
		String output = "\n";

		// output the pn(x)'s formula
		output += "P" + n + "(x)=f(x0)+";
		for (int i = 1; i < n + 1; i++) {
			for (int j = 0; j < i; j++)
				output += "(x-x" + j + ")";
			output += "f[";
			for (int j = 0; j < i + 1; j++) {
				output += "x(i+" + j + ")";
				if (j != i)
					output += ",";
			}
			output += "]";
			if (i == n)
				output += "=";
			else
				output += "+";
		}

		// output the formula that has values
		output += "\n" + "P" + n + "(" + x + ")=(" + y[0] + ")+";
		for (int i = 1; i < n + 1; i++) {
			for (int j = 0; j < i; j++)
				output += "(" + x + "-(" + this.x[j] + "))";
			output += "*(" + fmt.format(differenceQuotient(i, this.x, this.y))
					+ ")";
			if (i == n)
				output += "=";
			else
				output += "+";
		}

		// compute the result and output it
		double sum = 0;
		for (int i = 0; i < n + 1; i++) {
			double product = 1;
			for (int j = 0; j < i; j++)
				product *= x - this.x[j];
			sum += product * differenceQuotient(i, this.x, this.y);
		}
		output += "\n" + fmt.format(sum);

		output += "\n\n结果: " + fmt.format(sum);

		return output;
	}

	/**
	 * 计算一个区间的差分.定义n,赋值n=节点数减1,定义数组x[n+1]用于保存x[i]=xi(i=0,1,...,n);数组y同理.
	 * 注意:节点的xi应该等距;否则计算此区间的差分将无意义
	 * 
	 * @param order
	 *            差分的阶数
	 * @param subscript
	 *            差分的下标
	 * @param n
	 *            节点数减1的值
	 * @param x
	 *            x[i]=xi(i=0,1,...,n)
	 * @param y
	 *            y[i]=yi(i=0,1,...,n)
	 * @return 区间[x0,xn]的差分(阶数为order,下标为subscript)
	 */
	public static double difference(int order, int subscript, int n,
			double[] x, double[] y) {
		double result;

		if (order == 1)
			result = y[subscript + 1] - y[subscript];
		else
			result = difference(order - 1, subscript + 1, n, x, y)
					- difference(order - 1, subscript, n, x, y);

		return result;

	}

	/**
	 * 输出差分表.注意:节点的xi应该等距,否则就返回错误信息
	 * 
	 * @return 若节点的xi等距,则返回差分表;否则返回错误信息
	 */
	public String differenceChart() {
		String output = "\n";

		if (n > 1 && Math.abs((x[2] - x[1]) - (x[1] - x[0])) > 1e-6)
			output += "节点的xi不等间距!";
		else {
			output += "x\ty\t";
			for (int i = 1; i < n + 1; i++)
				output += "delta" + i + "y" + "\t";
			output += "\n";

			for (int i = 0; i < n + 1; i++) {
				output += x[i] + "\t" + y[i] + "\t";
				for (int j = 1; j < n + 1 - i; j++)
					output += fmt.format(difference(j, i, n, x, y)) + "\t";
				output += "\n";

			}
		}

		return output;
	}

	/**
	 * 牛顿前插公式,注意:节点的xi应该等距,否则就返回错误信息
	 * 
	 * @param x
	 *            Pn(x)的x
	 * @return 若节点xi等距则返回结果;否则返回错误信息
	 */
	public String newtonSequence(double x) {
		String output = "\n";

		if (n > 1
				&& Math.abs((this.x[2] - this.x[1]) - (this.x[1] - this.x[0])) > 1e-6)
			output += "节点的xi不等间距!";
		else {
			if (n == 0)
				output += "结果:" + this.y[n];
			else {
				double h = this.x[1] - this.x[0];
				double t = (x - this.x[0]) / h;

				// output t
				output += "t=(x-x0)/h=" + "((" + x + ")-(" + this.x[0] + "))/"
						+ h + "=" + t;
				output += "\n";

				// output the formula
				output += "\n" + "P" + n + "(x)=y0+";
				for (int i = 1; i < n + 1; i++) {
					output += "t";
					for (int j = 1; j < i; j++)
						output += "(t-" + j + ")";
					output += "*delta(" + i + ")yo/" + i + "!";
					if (i == n)
						output += "=";
					else
						output += "+";
				}

				// output the formula that has values
				output += "\n" + "P" + n + "(" + x + ")=(" + y[0] + ")+";
				for (int i = 1; i < n + 1; i++) {
					output += "(" + fmt.format(t) + ")";
					for (int j = 1; j < i; j++)
						output += "(" + fmt.format(t) + "-" + j + ")";
					output += "*"
							+ fmt.format(difference(i, 0, n, this.x, this.y))
							+ "/" + factorial(i);
					if (i == n)
						output += "=";
					else
						output += "+";
				}

				// compute and output the result
				double sum = this.y[0];
				for (int i = 1; i < n + 1; i++) {
					double product = t;
					for (int j = 1; j < i; j++)
						product *= t - j;
					sum += product * difference(i, 0, n, this.x, this.y)
							/ factorial(i);
				}
				output += "\n" + fmt.format(sum);

				output += "\n\n结果: " + fmt.format(sum);
			}
		}

		return output;
	}

	/**
	 * 牛顿后插公式,注意:节点的xi应该等距,否则就返回错误信息
	 * 
	 * @param x
	 *            Pn(x)的x
	 * @return 若节点xi等距则返回结果;否则返回错误信息
	 */
	public String newtonInverse(double x) {
		String output = "\n";

		if (n > 1
				&& Math.abs((this.x[2] - this.x[1]) - (this.x[1] - this.x[0])) > 1e-6)
			output += "节点的xi不等间距!";
		else {
			if (n == 0)
				output += "结果:" + this.y[n];
			else {
				double h = this.x[1] - this.x[0];
				double t = (x - this.x[n]) / h;

				// output t
				output += "t=(x-xn)/h=" + "((" + x + ")-(" + this.x[n] + "))/"
						+ h + "=" + t;
				output += "\n";

				// output the formula
				output += "\n" + "P" + n + "(x)=yn+";
				for (int i = 1; i < n + 1; i++) {
					output += "t";
					for (int j = 1; j < i; j++)
						output += "(t+" + j + ")";
					output += "*delta(" + i + ")y(" + n + "-" + i + ")/" + i
							+ "!";
					if (i == n)
						output += "=";
					else
						output += "+";
				}

				// output the formula that has values
				output += "\n" + "P" + n + "(" + x + ")=(" + y[n] + ")+";
				for (int i = 1; i < n + 1; i++) {
					output += "(" + fmt.format(t) + ")";
					for (int j = 1; j < i; j++)
						output += "(" + fmt.format(t) + "+" + j + ")";
					output += "*"
							+ fmt
									.format(difference(i, n - i, n, this.x,
											this.y)) + "/" + factorial(i);
					if (i == n)
						output += "=";
					else
						output += "+";
				}

				// compute and output the result
				double sum = this.y[n];
				for (int i = 1; i < n + 1; i++) {
					double product = t;
					for (int j = 1; j < i; j++)
						product *= t + j;
					sum += product * difference(i, n - i, n, this.x, this.y)
							/ factorial(i);
				}
				output += "\n" + fmt.format(sum);

				output += "\n\nso the result is " + fmt.format(sum);
			}
		}

		return output;
	}

	/**
	 * 拉格朗日插值公式
	 * 
	 * @param x
	 *            Ln(x)的x
	 * @return 结果
	 */
	public String lagrange(double x) {
		String output = "\n";

		// output the formula
		output += "L" + n + "(x)=";
		for (int i = 0; i < n + 1; i++) {
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					output += "(x-x" + j + ")";
			output += "*f(x" + i + ")/";
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					output += "(x" + i + "-x" + j + ")";
			if (i == n)
				output += "=";
			else
				output += "+";
		}
		output += "\n";

		// output the formula that has values
		output += "L" + n + "(" + x + ")=";
		for (int i = 0; i < n + 1; i++) {
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					output += "(" + x + "-" + this.x[j] + ")";
			output += "*(" + this.y[i] + ")/";
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					output += "(" + this.x[i] + "-" + this.x[j] + ")";
			if (i == n)
				output += "=";
			else
				output += "+";
		}
		output += "\n";

		// compute and output the result
		double sum = 0;
		for (int i = 0; i < n + 1; i++) {
			double numerator = 1;
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					numerator *= x - this.x[j];
			numerator *= this.y[i];

			double denominator = 1;
			for (int j = 0; j < n + 1; j++)
				if (j != i)
					denominator *= this.x[i] - this.x[j];

			sum += numerator / denominator;
		}
		output += fmt.format(sum);

		output += "\n\n结果: " + fmt.format(sum);

		return output;
	}

	/**
	 * 等距节点的正函数反插值(书中公式(5.104)),注意:只有x0,x1,...,xn等距才适用,否则返回错误信息
	 * 
	 * @param y
	 *            y=Pn(x)的y
	 * @return 若节点xi等距则返回结果;否则返回错误信息
	 */
	public String sameDistanceInverseIntepolation(double y) {
		String output = "\n";

		double t, tLast;

		if (n == 0)
			output += "结果:" + this.x[n];
		else {
			if (n > 1 && Math.abs((x[2] - x[1]) - (x[1] - x[0])) > 1e-6)
				output += "节点的xi不等间距!";
			else {
				// output the formula with value
				output += "h=" + fmt.format(Math.abs((x[1] - x[0])));
				output += "\nt=(" + y + "-" + this.y[0] + ")/"
						+ fmt.format(difference(1, 0, n, this.x, this.y)) + "-";
				for (int i = 2; i < n + 1; i++) {
					output += fmt.format(difference(i, 0, n, this.x, this.y))
							+ "t";
					for (int j = 1; j < i; j++)
						output += "(t-" + j + ")";
					output += "/" + factorial(i) + "*"
							+ fmt.format(difference(1, 0, n, this.x, this.y));
					if (i != n)
						output += "-";
				}

				t = (y - this.y[0]) / difference(1, 0, n, this.x, this.y);
				output += "\nt0=" + t + "\n";
				int counter = 1;
				do {
					tLast = t;
					t = (y - this.y[0]) / difference(1, 0, n, this.x, this.y);
					for (int i = 2; i < n + 1; i++) {
						double product = 1, quotient = 1;
						product *= difference(i, 0, n, this.x, this.y) * tLast;
						for (int j = 1; j < i; j++)
							product *= tLast - j;
						quotient += factorial(i)
								* difference(1, 0, n, this.x, this.y);
						t -= product / quotient;
					}
					output += "t" + (counter++) + "=" + fmt.format(t) + "\n";
				} while (Math.abs(t - tLast) > 1e-5);

				output += "x=x0+th=" + x[0] + "+" + t + "*"
						+ fmt.format(Math.abs((x[1] - x[0]))) + "="
						+ fmt.format(x[0] + t * Math.abs((x[1] - x[0])));
			}
		}

		return output;

	}

	/**
	 * 计算阶乘
	 * 
	 * @param n
	 *            计算n的阶乘,n必须为不小于0的整数
	 * @return
	 */
	public static int factorial(int n) {
		int result;

		if (n == 0)
			result = 1;
		else
			result = n * factorial(n - 1);

		return result;
	}

}