integration.java

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

				double i = (rightEndpoint - leftEndpoint)
						* (fx[0] / 8 + 3 * fx[1] / 8 + 3 * fx[2] / 8 + fx[3] / 8);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 4: {
				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(7/90*" + fmt.format(fx[0]) + "+";
				output += "\n\n16/45*" + fmt.format(fx[1]) + "+";
				output += "\n\n2/15*" + fmt.format(fx[2]) + "+";
				output += "\n\n16/45*" + fmt.format(fx[3]) + "+";
				output += "\n\n7/90*" + fmt.format(fx[4]) + ")=";
				double i = (rightEndpoint - leftEndpoint)
						* (7 * fx[0] / 90 + 16 * fx[1] / 45 + 2 * fx[2] / 15
								+ 16 * fx[3] / 45 + 7 * fx[4] / 90);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 5: {
				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(19/288*" + fmt.format(fx[0]) + "+";
				output += "\n\n25/96*" + fmt.format(fx[1]) + "+";
				output += "\n\n25/144*" + fmt.format(fx[2]) + "+";
				output += "\n\n25/144*" + fmt.format(fx[3]) + "+";
				output += "\n\n25/96*" + fmt.format(fx[4]) + "+";
				output += "\n\n19/288*" + fmt.format(fx[5]) + ")=";
				double 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);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 6: {
				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(41/840*" + fmt.format(fx[0]) + "+";
				output += "\n\n9/35*" + fmt.format(fx[1]) + "+";
				output += "\n\n9/280*" + fmt.format(fx[2]) + "+";
				output += "\n\n34/105*" + fmt.format(fx[3]) + "+";
				output += "\n\n9/280*" + fmt.format(fx[4]) + "+";
				output += "\n\n9/35*" + fmt.format(fx[5]) + "+";
				output += "\n\n41/840*" + fmt.format(fx[6]) + ")=";
				double 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);
				output += "\n\n" + fmt.format(i);
				break;
			}
			case 7: {
				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(751/17280*" + fmt.format(fx[0]) + "+";
				output += "\n\n3577/17280*" + fmt.format(fx[1]) + "+";
				output += "\n\n1323/17280*" + fmt.format(fx[2]) + "+";
				output += "\n\n2989/17280*" + fmt.format(fx[3]) + "+";
				output += "\n\n2989/17280*" + fmt.format(fx[4]) + "+";
				output += "\n\n1323/17280*" + fmt.format(fx[5]) + "+";
				output += "\n\n3577/17280*" + fmt.format(fx[6]) + "+";
				output += "\n\n751/17280*" + fmt.format(fx[7]) + ")=";
				double 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);
				output += "\n\n" + fmt.format(i);
				break;
			}
			}
		} else
			output += "n>7,不建议用牛顿柯特斯公式求积分.";

		return output;
	}

	/**
	 * 复化牛顿柯特斯公式:积分区间分为m大段,每大段都分为n小段并用牛顿柯特斯公式.
	 * 
	 * @param m
	 *            大段数目
	 * @param n
	 *            每大段分为n小段,用牛顿柯特斯公式
	 * @return 积分值
	 */
	public double minuteNewtonCotes(int m, int n) {
		double h;// 大段步长
		double[] x;// 把积分区间分为m大段的节点
		Integration[] problems;// 每大段都构造一个新积分,然后用牛顿柯特斯公式
		double I=0;//积分值
		
		h = (rightEndpoint - leftEndpoint) / m;
		x = new double[m + 1];
		for (int i = 0; i < m + 1; i++)
			x[i] = leftEndpoint + i * h;// (i=0,1,...,m)

		problems = new Integration[m];

		for (int i = 0; i < m; i++)
			problems[i] = new Integration(x[i], x[i + 1], function);// 构造每大段的积分

		for (int i = 0; i < m; i++)
			I+=problems[i].getSoutionByNewtonCotes(n);// 每大段用牛顿柯特斯公式
		
		return I;
	}

	/**
	 * 龙贝格法,用本方法求积分非常好用.
	 * 
	 * @param errorAllowed
	 *            允许误差
	 * @return 将会输出用龙贝格法的计算过程和积分结果I
	 */
	public String romberg(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;
		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)

		int i = 0;
		boolean stable = false;// 是否稳定至要求精度
		while (!stable) {
			// 计算T2^i
			if (i == 0) {
				output += "\n\nT" + (int) Math.pow(2, i) + "=0.5*[f(" + x[0]
						+ ")+f(" + x[n] + ")]=";
			} else {
				output += "\n\nT" + (int) Math.pow(2, i) + "=0.5*T"
						+ (int) Math.pow(2, i - 1) + "+(" + rightEndpoint + "-"
						+ leftEndpoint + ")/" + (int) Math.pow(2, i) + "*[";
				for (int k = 1; k <= (int) Math.pow(2, i - 1); k++) {
					output += "f(" + (2 * k - 1) + "/" + (int) Math.pow(2, i)
							+ ")";
					if (k != (int) Math.pow(2, i - 1))
						output += "+";
				}
				output += "]=";
			}
			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;
			}
			output += fmt.format(t[i]);

			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)
				output += "\n\nS" + (int) Math.pow(2, i - 1) + "=4*T"
						+ (int) Math.pow(2, i) + "/3" + "-T"
						+ (int) Math.pow(2, i - 1) + "/3=";
				s[i - 1] = 4 * t[i] / 3 - t[i - 1] / 3;
				output += fmt.format(s[i - 1]);

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

			if (!stable && i > 1) {// 若是第二行则只有T2,S2,否则若还不稳定则计算C2^(i-2)
				output += "\n\nC" + (int) Math.pow(2, i - 2) + "=4^2*S"
						+ (int) Math.pow(2, i - 1) + "/(4^2-1)" + "-S"
						+ (int) Math.pow(2, i - 2) + "/(4^2-1)=";
				c[i - 2] = 16 * s[i - 1] / 15 - s[i - 2] / 15;
				output += fmt.format(c[i - 2]);

				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)
				output += "\n\nR" + (int) Math.pow(2, i - 3) + "=4^3*C"
						+ (int) Math.pow(2, i - 2) + "/(4^3-1)" + "-C"
						+ (int) Math.pow(2, i - 3) + "/(4^3-1)=";
				r[i - 3] = 64 * c[i - 2] / 63 - c[i - 3] / 63;
				output += fmt.format(r[i - 3]);

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

			}

			i++;
		}

		output += "\n\nI=" + fmt.format(I);

		return output;
	}
}