integration.java

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

package numericalAnalysis.integration;

import java.text.DecimalFormat;
import numericalAnalysis.function.Function;

/**
 * 积分求解,包括牛顿柯特斯公式和龙贝格法.对照书:数值分析第二版,史万明等编,北京理工大学出版社
 * 
 * @author 山
 * 
 */
public class Integration {
	private double leftEndpoint, rightEndpoint;// 积分区间的左端点,右端点
	private Function function;// 被积函数所属类
	private DecimalFormat fmt = new DecimalFormat("0.######");// 输出六位小数

	/**
	 * 新建一个积分求解问题.外部类中必须定义一个类实现Function接口,这个类的对象的f方法即被积函数f(x).
	 * 
	 * @param leftEndpoint
	 *            积分区间的左端点
	 * @param rightEndpoint
	 *            积分区间的右端点
	 * @param function
	 *            function的f方法即为被积函数
	 */
	public Integration(double leftEndpoint, double rightEndpoint,
			Function function) {
		this.leftEndpoint = leftEndpoint;
		this.rightEndpoint = rightEndpoint;
		this.function = function;
	}

	/**
	 * 获取积分结果(用龙贝格法求积分)
	 * 
	 * @param errorAllowed
	 *            允许误差
	 * @return 积分结果
	 */
	public double getSolution(double errorAllowed) {
		final int n = 8192;// 这里为了使龙贝格法达到要求允许误差,建议分段数n为8192(即2^13)(利用循环容易实现x与fx的赋值,无论n多大)
		double h;// 步长
		double[] x, fx;// 节点,节点代入被积函数的值
		double[] t = new double[20], s = new double[20], c = new double[20], r = new double[20];// 书上的T,S,C,R
		double I = 0;

		h = (rightEndpoint - leftEndpoint) / n;
		x = new double[n + 1];
		fx = new double[n + 1];
		for (int i = 0; i < n + 1; i++)
			x[i] = leftEndpoint + i * h;// (i=0,1,...,n)
		for (int i = 0; i < n + 1; i++)
			fx[i] = function.f(x[i]);// (i=0,1,...,n)

		int i = 0;
		boolean stable = false;// 是否稳定至要求精度
		while (!stable) {
			// 计算T2^i
			double sum = 0;
			if (i == 0)
				t[i] = (rightEndpoint - leftEndpoint) / 2 * (fx[0] + fx[n]);
			else {
				for (int k = 1; k <= (int) Math.pow(2, i - 1); k++)
					sum += fx[(int) (n / (int) Math.pow(2, i) * (2 * k - 1))];
				t[i] = 0.5 * t[i - 1] + (rightEndpoint - leftEndpoint)
						/ (int) Math.pow(2, i) * sum;
			}

			if (i == 1) {
				if (Math.abs(t[i] - t[i - 1]) < errorAllowed) {
					I = t[i];
					stable = true;
				}
			} else if (i == 2) {
				if (Math.abs(t[i] - s[0]) < errorAllowed) {
					I = t[i];
					stable = true;
				}
			} else if (i == 3) {
				if (Math.abs(t[i] - c[0]) < errorAllowed) {
					I = t[i];
					stable = true;
				}
			} else if (i > 3) {
				if (Math.abs(t[i] - r[i - 4]) < errorAllowed) {
					I = t[i];
					stable = true;
				}
			}

			if (i > 0 && !stable) {// 若是第一行则只有T1,否则计算S2^(i-1)
				s[i - 1] = 4 * t[i] / 3 - t[i - 1] / 3;

				if (Math.abs(t[i] - s[i - 1]) < errorAllowed) {
					I = s[i - 1];
					stable = true;
				}// 已经稳定
			}

			if (!stable && i > 1) {// 若是第二行则只有T2,S2,否则若还不稳定则计算C2^(i-2)
				c[i - 2] = 16 * s[i - 1] / 15 - s[i - 2] / 15;

				if (Math.abs(s[i - 1] - c[i - 2]) < errorAllowed) {
					I = c[i - 2];
					stable = true;
				}// 已经稳定
			}

			if (!stable && i > 2) {// 若是第三行则只有T2,S2,C1,否则若还不稳定则计算R2^(i-3)
				r[i - 3] = 64 * c[i - 2] / 63 - c[i - 3] / 63;

				if (Math.abs(c[i - 2] - r[i - 3]) < errorAllowed) {
					I = r[i - 3];
					stable = true;
				}// 稳定
			}

			i++;
		}

		return I;
	}

	/**
	 * 求解积分(用牛顿柯特斯公式)
	 * 
	 * @param n 积分区间分为n段(n<=7)
	 * @return 积分值
	 */
	public double getSoutionByNewtonCotes(int n) {
		double h;// 步长
		double[] x, fx;// 节点,节点代入被积函数的值
		double I = 0;// 积分值

		h = (rightEndpoint - leftEndpoint) / n;
		x = new double[n + 1];
		fx = new double[n + 1];
		for (int i = 0; i < n + 1; i++)
			x[i] = leftEndpoint + i * h;// (i=0,1,...,n)
		for (int i = 0; i < n + 1; i++)
			fx[i] = function.f(x[i]);// (i=0,1,...,n)

		switch (n) {
		case 1: {
			I = (rightEndpoint - leftEndpoint) * (fx[0] / 2 + fx[1] / 2);
			break;
		}
		case 2: {
			I = (rightEndpoint - leftEndpoint)
					* (fx[0] / 6 + 4 * fx[1] / 6 + fx[2] / 6);
			break;
		}
		case 3: {
			I = (rightEndpoint - leftEndpoint)
					* (fx[0] / 8 + 3 * fx[1] / 8 + 3 * fx[2] / 8 + fx[3] / 8);
			break;
		}
		case 4: {
			I = (rightEndpoint - leftEndpoint)
					* (7 * fx[0] / 90 + 16 * fx[1] / 45 + 2 * fx[2] / 15 + 16
							* fx[3] / 45 + 7 * fx[4] / 90);
			break;
		}
		case 5: {
			I = (rightEndpoint - leftEndpoint)
					* (19 * fx[0] / 288 + 25 * fx[1] / 96 + 25 * fx[2] / 144
							+ 25 * fx[3] / 144 + 25 * fx[4] / 96 + 19 * fx[5] / 288);
			break;
		}
		case 6: {
			I = (rightEndpoint - leftEndpoint)
					* (41 * fx[0] / 840 + 9 * fx[1] / 35 + 9 * fx[2] / 280 + 34
							* fx[3] / 105 + 9 * fx[4] / 280 + 9 * fx[5] / 35 + 41 * fx[6] / 840);
			break;
		}
		case 7: {
			I = (rightEndpoint - leftEndpoint)
					* (751 * fx[0] / 17280 + 3577 * fx[1] / 17280 + 1323
							* fx[2] / 17280 + 2989 * fx[3] / 17280 + 2989
							* fx[4] / 17280 + 1323 * fx[5] / 17280 + 3577
							* fx[6] / 17280 + 751 * fx[6] / 17280);
			break;
		}
		}

		return I;
	}

	/**
	 * 牛顿柯特斯公式(n<=7)
	 * 
	 * @param n
	 *            积分区间分为n段(n<=7)
	 * @return 将会输出步长h,各节点xi,各节点代入被积函数的值f(xi),积分I
	 */
	public String newtonCotes(int n) {
		double h;// 步长
		double[] x, fx;// 节点,节点代入被积函数的值
		String output = "\n";

		h = (rightEndpoint - leftEndpoint) / n;
		x = new double[n + 1];
		fx = new double[n + 1];
		for (int i = 0; i < n + 1; i++)
			x[i] = leftEndpoint + i * h;// (i=0,1,...,n)
		for (int i = 0; i < n + 1; i++)
			fx[i] = function.f(x[i]);// (i=0,1,...,n)

		if (n <= 7) {
			switch (n) {
			case 1: {
				output += "h=" + fmt.format(h);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "x" + i + "=" + fmt.format(x[i]);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "f(x" + i + ")=" + fx[i];
				output += "\n\nI=(" + fmt.format(rightEndpoint) + "-"
						+ fmt.format(leftEndpoint) + ")*";
				output += "\n\n(1/2*" + fmt.format(fx[0]) + "+";
				output += "\n\n1/2*" + fmt.format(fx[1]) + ")=";
				double i;
				i = (rightEndpoint - leftEndpoint) * (fx[0] / 2 + fx[1] / 2);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 2: {
				output += "h=" + fmt.format(h);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "x" + i + "=" + fmt.format(x[i]);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "f(x" + i + ")=" + fmt.format(fx[i]);
				output += "\n\nI=(" + fmt.format(rightEndpoint) + "-"
						+ fmt.format(leftEndpoint) + ")*";
				output += "\n\n(1/6*" + fmt.format(fx[0]) + "+";
				output += "\n\n4/6*" + fmt.format(fx[1]) + "+";
				output += "\n\n1/6*" + fmt.format(fx[2]) + ")=";
				double i = (rightEndpoint - leftEndpoint)
						* (fx[0] / 6 + 4 * fx[1] / 6 + fx[2] / 6);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 3: {
				output += "h=" + fmt.format(h);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "x" + i + "=" + fmt.format(x[i]);
				for (int i = 0; i < n + 1; i++)
					output += "\n\n" + "f(x" + i + ")=" + fmt.format(fx[i]);
				output += "\n\nI=(" + fmt.format(rightEndpoint) + "-"
						+ fmt.format(leftEndpoint) + ")*";
				output += "\n\n(1/8*" + fmt.format(fx[0]) + "+";
				output += "\n\n3/8*" + fmt.format(fx[1]) + "+";
				output += "\n\n3/8*" + fmt.format(fx[2]) + "+";
				output += "\n\n1/8*" + fmt.format(fx[3]) + ")=";