linearequationgroup.java

数学：数值分析：线性方程组 线性方程组的多种解法

					xLast[i] = x[i];
				for (int i = 1; i < n + 1; i++) {
					double sum = 0;
					for (int j = 1; j < n + 1; j++)
						sum += g[i][j] * x[j];
					x[i] = sum + f[i];
				}
				double[] xsubtract = new double[n + 1];
				for (int i = 1; i < n + 1; i++)
					xsubtract[i] = x[i] - xLast[i];
				xNorm = Norm.vectorOne(n, xsubtract);
			} while (xNorm > (1 - gNorm) * errAllowed / gNorm);
			return output;
		} else {
			output += "\nusing Jacobi iterative method is not convergent.\n";
			return output;
		}
	}

	// ----------------------------------------------------------------------------------------
	// Gauss-Seidel iterative method(input the error allowed and the initialized
	// value).
	// ----------------------------------------------------------------------------------------
	public String gaussSeidelIterative(double errAllowed, double[] xIni) {
		double[][] m1 = new double[n + 1][n + 1], m2 = new double[n + 1][n + 1];
		double[] nVector = new double[n + 1], x = new double[n + 1], xLast = new double[n + 1];
		String output = "\n";
		boolean convergent = true;

		// judge whether using this method is convergent(not right but useful)
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				if (a[i][j] == 0)
					convergent = false;
		for (int i = 1; i < n + 1; i++) {
			double sum = 0;
			for (int j = 1; j < n + 1; j++)
				if (j != i)
					sum += Math.abs(a[i][j]);
			if (Math.abs(a[i][i]) <= sum)
				convergent = false;
		}

		if (convergent == true) {// compute the M1(-D-1L) matrix,the
			// M2(-D-1U) matrix and the N(D-1B)
			// vector then output them
			output += "\n\nthe M1 matrix is:\n";
			for (int i = 1; i < n + 1; i++)
				for (int j = 1; j < n + 1; j++)
					if (i > j) {
						m1[i][j] = (-1) * a[i][j] / a[i][i];
						output += "m" + i + j + "=(" + fmt.format(m1[i][j])
								+ ")   ";
					} else {
						m1[i][j] = 0;
						output += "m" + i + j + "=(" + fmt.format(m1[i][j])
								+ ")   ";
						if (j == n)
							output += "\n";
					}
			output += "\n\nthe M2 matrix is:\n";
			for (int i = 1; i < n + 1; i++)
				for (int j = 1; j < n + 1; j++)
					if (i < j) {
						m2[i][j] = (-1) * a[i][j] / a[i][i];
						output += "m" + i + j + "=(" + fmt.format(m2[i][j])
								+ ")   ";
						if (j == n)
							output += "\n";
					} else {
						m2[i][j] = 0;
						output += "m" + i + j + "=(" + fmt.format(m2[i][j])
								+ ")   ";
					}
			output += "\n\nthe N vector is:\n";
			for (int i = 1; i < n + 1; i++) {
				nVector[i] = b[i] / a[i][i];
				output += "n" + i + "=(" + fmt.format(nVector[i]) + ")";
			}

			// output the iterative formula
			output += "\n\nthe iterative formula is:\n";
			for (int i = 1; i < n + 1; i++)
				for (int j = 0; j < n + 2; j++)
					if (j == 0)
						output += "x" + i + "(k+1)=";
					else if (j == n + 1)
						output += "(" + fmt.format(nVector[i]) + ")\n";
					else if (i < j)
						output += "(" + fmt.format(m2[i][j]) + "x" + j
								+ "(k))+";
					else
						output += "(" + fmt.format(m1[i][j]) + "x" + j
								+ "(k+1))+";

			// compute xi and output xi
			output += "\n\nthe result is:\n";
			double xNorm = 0;
			int count = 0;
			for (int i = 1; i < n + 1; i++)
				xLast[i] = x[i] = xIni[i];
			do {
				for (int i = 1; i < n + 1; i++)
					output += "x" + i + "(" + count + ")=(" + fmt.format(x[i])
							+ ")\n";
				count++;
				output += "\n";
				for (int i = 1; i < n + 1; i++)
					xLast[i] = x[i];
				for (int i = 1; i < n + 1; i++) {
					double sum1 = 0, sum2 = 0;
					for (int j = 1; j < i; j++)
						sum1 += m1[i][j] * x[j];
					for (int j = i + 1; j < n + 1; j++)
						sum2 += m2[i][j] * xLast[j];
					x[i] = sum1 + sum2 + nVector[i];
				}
				double[] xsubtract = new double[n + 1];
				for (int i = 1; i < n + 1; i++)
					xsubtract[i] = x[i] - xLast[i];
				xNorm = Norm.vectorOne(n, xsubtract);
			} while (xNorm > errAllowed);
			return output;
		} else {
			output += "\nusing Gauss-Seidel iterative method is not convergent.\n";
			return output;
		}
	}

	// ----------------------------------------------------------------------------------------
	// Relaxation iterative method(input the error allowed,the initialized
	// value and the relaxation parameter).
	// ----------------------------------------------------------------------------------------
	public String relaxationIterative(double errAllowed, double[] xIni,
			double omega) {
		double[][] m1 = new double[n + 1][n + 1], m2 = new double[n + 1][n + 1];
		double[] nVector = new double[n + 1], x = new double[n + 1], xLast = new double[n + 1];
		String output = "\n";

		// compute the M1(-D-1L) matrix,the M2(-D-1U) matrix and the N(D-1B)
		// vector then output them
		output += "\n\nthe M1 matrix is:\n";
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				if (i > j) {
					m1[i][j] = omega * (-1) * a[i][j] / a[i][i];
					output += "m" + i + j + "=(" + fmt.format(m1[i][j])
							+ ")   ";
				} else {
					m1[i][j] = 0;
					output += "m" + i + j + "=(" + fmt.format(m1[i][j])
							+ ")   ";
					if (j == n)
						output += "\n";
				}
		output += "\n\nthe M2 matrix is:\n";
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				if (i < j) {
					m2[i][j] = omega * (-1) * a[i][j] / a[i][i];
					output += "m" + i + j + "=(" + fmt.format(m2[i][j])
							+ ")   ";
					if (j == n)
						output += "\n";
				} else {
					m2[i][j] = 0;
					output += "m" + i + j + "=(" + fmt.format(m2[i][j])
							+ ")   ";
				}
		output += "\n\nthe N vector is:\n";
		for (int i = 1; i < n + 1; i++) {
			nVector[i] = omega * b[i] / a[i][i];
			output += "n" + i + "=(" + fmt.format(nVector[i]) + ")";
		}

		// output the iterative formula
		output += "\n\nthe iterative formula is:\n";
		for (int i = 1; i < n + 1; i++)
			for (int j = 0; j < n + 2; j++)
				if (j == 0)
					output += "x" + i + "(k+1)=";
				else if (j == n + 1)
					output += "(" + fmt.format(nVector[i]) + ")\n";
				else if (i < j)
					output += "(" + fmt.format(m2[i][j]) + "x" + j + "(k))+";
				else if (i == j)
					output += "(" + fmt.format(1 - omega) + "x" + j + "(k))+";
				else
					output += "(" + fmt.format(m1[i][j]) + "x" + j + "(k+1))+";

		// compute xi and output xi
		output += "\n\nthe result is:\n";
		double xNorm = 0;
		int count = 0;
		for (int i = 1; i < n + 1; i++)
			xLast[i] = x[i] = xIni[i];
		do {
			for (int i = 1; i < n + 1; i++)
				output += "x" + i + "(" + count + ")=(" + fmt.format(x[i])
						+ ")\n";
			count++;
			output += "\n";
			for (int i = 1; i < n + 1; i++)
				xLast[i] = x[i];
			for (int i = 1; i < n + 1; i++) {
				double sum1 = 0, sum2 = 0;
				for (int j = 1; j < i; j++)
					sum1 += m1[i][j] * x[j];
				for (int j = i + 1; j < n + 1; j++)
					sum2 += m2[i][j] * xLast[j];
				x[i] = sum1 + sum2 + nVector[i] + (1 - omega) * xLast[i];
			}
			double[] xsubtract = new double[n + 1];
			for (int i = 1; i < n + 1; i++)
				xsubtract[i] = x[i] - xLast[i];
			xNorm = Norm.vectorInfinitive(n, xsubtract);
		} while (xNorm > errAllowed);

		return output;
	}

	// ----------------------------------------------------------------------------------------
	// Main element method.
	// ----------------------------------------------------------------------------------------
	public String mainElement() {
		double[][] aElim = new double[n + 1][n + 1];
		double[] bElim = new double[n + 1], x = new double[n + 1];
		String output = "\n";

		// initialize the aElim matrix and bElim vector
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				aElim[i][j] = a[i][j];
		for (int i = 1; i < n + 1; i++)
			bElim[i] = b[i];

		// elimination and output the equation group after elimination
		output += "\n\nthe equation group after elimination:\n";
		for (int i = 1; i < n + 1; i++) {
			double mainElement;
			int mainElementMark;

			// select the main element of the line i,stored in row maxMark
			mainElement = aElim[i][i];
			mainElementMark = i;
			for (int j = i + 1; j < n + 1; j++)
				if (Math.abs(aElim[j][i]) > Math.abs(mainElement)) {
					mainElement = aElim[j][i];
					mainElementMark = j;
				}

			// swap all the elements of row i and row mainElementMark
			for (int j = i; j < n + 2; j++)
				if (j == n + 1) {
					double temp;
					temp = bElim[mainElementMark];
					bElim[mainElementMark] = bElim[i];
					bElim[i] = temp;
				} else {
					double temp;
					temp = aElim[mainElementMark][j];
					aElim[mainElementMark][j] = aElim[i][j];
					aElim[i][j] = temp;
				}

			// output the equation i after selecting the main element
			for (int j = i; j < n + 2; j++)
				if (j == n + 1)
					output += "=(" + bElim[i] + ")\n";
				else {
					output += "(" + fmt.format(aElim[i][j]) + "x" + j + ")";
					if (j != n)
						output += "+";
				}

			// elimination of the matrix whose level is i
			double l;
			for (int j = i + 1; j < n + 1; j++) {
				l = aElim[j][i] / aElim[i][i];
				for (int k = i; k < n + 2; k++)
					if (k == n + 1)
						bElim[j] = bElim[j] - l * bElim[i];
					else
						aElim[j][k] = aElim[j][k] - l * aElim[i][k];
			}
		}

		// compute and output xi
		output += "\n\nthe result is:\n";
		for (int i = n; i > 0; i--) {
			double sum = 0;
			for (int k = i + 1; k < n + 1; k++)
				sum += aElim[i][k] * x[k];
			x[i] = (bElim[i] - sum) / aElim[i][i];
			output += "x" + i + "=(" + x[i] + ")\n";
		}

		return output;
	}
}