differential.java

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

package numericalAnalysis.differential;

import java.text.DecimalFormat;

/**
 * 微分求解,包括牛顿柯特斯公式和龙贝格法.本类不容易理解,使用本类前必读懂书上197页例6.12!对照书:数值分析第二版,史万明等编,北京理工大学出版社.
 * 
 * @author 山
 * 
 */
public class Differential {
	private static DecimalFormat fmt = new DecimalFormat("0.######");

	/**
	 * 利用n点式求解微分f'(x),f''(x),...直至f(k)(x).参考例题:书上P197,例6.12
	 * 
	 * @param order
	 *            求解微分f'(x),f''(x),...直至f(k)(x)的最高阶数k
	 * @param n
	 *            n点式的n,如:表格函数中有4个点,则为x0,x1,x2,x3,因此n=3.这(n+1)个点最后一个为f(k)(x)的x
	 * @param x
	 *            在外部类中定义x=new
	 *            double[n+1],x[i](i=0,1,...,n-1)为表格函数中除了x外的其他点,
	 *            由用户输入,x[n]必须为x
	 * @param y
	 *            y[i]为x[i]的对应y值
	 * @return 结果
	 */
	public static String differential(int order,int n, double[] x,
			double[] y) {
		String output = "\n";
		double[][] chart = new double[n+1+order][n+2];
		// 如书上197页的6.8表

		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; i++) {
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\t";
			chart[i][0] = x[i];
			chart[i][1] = y[i];
			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];
				}
				double getValue = differenceQuotient(j, tempX, tempY);
				chart[i][j + 1] = getValue;
				output += fmt.format(getValue) + "\t";
			}
			output += "\n";
		}

		for (int i = n+1; i < n+1+order; i++) {
			chart[i][0] = x[n];
			chart[i][1] = y[n];
			chart[i][n+1] = chart[n][n+1];
		}

		for (int j = n; j >= 2; j--)
			for (int i = n+1; i <= n + order; i++)
				chart[i][j] = chart[i - 1][j] + (chart[i][0] - chart[i - j][0])
						* chart[i][j + 1];

		for(int i=1;i<=order;i++){
		output += "\n\nf(" + i + ")(" + x[n] + ")=" + i + "!*P"
				+ n + "[";
		for (int j = 1; j <=i+1; j++) {
			output += x[n];
			if (j != i + 1)
				output += ",";
		}
		output += "]=";

		output += fmt.format(factorial(i)
				* chart[n+order][1+i]);
		}
		return output;
	}

	/**
	 * 计算一个区间的差商.定义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;
	}

	/**
	 * 计算阶乘
	 * 
	 * @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;
	}

}