linearequationgroup.java

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

			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 G matrix and the F vector then output them
			output += "\n\n矩阵G:\n";
			for (int i = 1; i < n + 1; i++)
				for (int j = 1; j < n + 1; j++) {
					if (i == j) {
						g[i][j] = 0;
						output += "g" + i + j + "=(" + fmt.format(g[i][j])
								+ ")   ";
					} else {
						g[i][j] = (-1) * a[i][j] / a[i][i];
						output += "g" + i + j + "=(" + fmt.format(g[i][j])
								+ ")   ";
						if (j == n)
							output += "\n";
					}
				}
			output += "\n\n向量F:\n";
			for (int i = 1; i < n + 1; i++) {
				f[i] = b[i] / a[i][i];
				output += "f" + i + "=(" + fmt.format(f[i]) + ")   ";
			}

			// output the iterative formula
			output += "\n\n迭代公式:\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(f[i]) + ")\n";
					else
						output += "(" + fmt.format(g[i][j]) + "x" + j + "(k))+";

			// compute xi and output xi
			output += "\n\n解:\n";
			double gNorm = Norm.matrixOne(n, g), 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 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 += "\n使用本方法不收敛\n";
			return output;
		}
	}

	/**
	 * 高斯-赛德尔迭代法(收敛判断未完成).
	 * 
	 * @param errAllowed
	 *            允许误差
	 * @param xIni
	 *            初值.注意:必须定义xIni[i]为xi.即xIni[0]不使用
	 * @return 结果
	 */
	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\n矩阵M1:\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\n矩阵M2:\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\n向量N:\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\n迭代公式:\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\n解:\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 += "\n使用本方法不收敛\n";
			return output;
		}
	}

	/**
	 * 松弛迭代法(收敛判断未完成).
	 * 
	 * @param errAllowed
	 *            允许误差
	 * @param xIni
	 *            初值.注意:必须定义xIni[i]为xi.即xIni[0]不使用
	 * @param factor
	 *            松弛因子
	 * @return 结果
	 */
	public String relaxationIterative(double errAllowed, double[] xIni,
			double factor) {
		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\n矩阵M1:\n";
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				if (i > j) {
					m1[i][j] = factor * (-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\n矩阵M2:\n";
		for (int i = 1; i < n + 1; i++)
			for (int j = 1; j < n + 1; j++)
				if (i < j) {
					m2[i][j] = factor * (-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\n向量N:\n";
		for (int i = 1; i < n + 1; i++) {
			nVector[i] = factor * b[i] / a[i][i];
			output += "n" + i + "=(" + fmt.format(nVector[i]) + ")";
		}

		// output the iterative formula
		output += "\n\n迭代公式:\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 - factor) + "x" + j + "(k))+";
				else
					output += "(" + fmt.format(m1[i][j]) + "x" + j + "(k+1))+";

		// compute xi and output xi
		output += "\n\n解:\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 - factor) * 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;
	}

	/**
	 * 列主元素消元法
	 * 
	 * @return 结果
	 */
	public String mainElementOfColumns() {
		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\n消元后的方程组:\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\n解:\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;
	}

	/**
	 * 求解过定方程组
	 * 
	 * @return 过定方程组的解
	 */
	public String overdetermined(){
		double[][] aT=new double[a[1].length][a.length],aM=new double[a[1].length][a[1].length];
		//a=A,aT=A',aM=A'A
		double[] ABM=new double[a[1].length];//b=B,ABM=A'B
		String output="\n";
		
		for(int i=1;i<aT.length;i++)
			for(int j=1;j<aT[1].length;j++)
				aT[i][j]=a[j][i];
		
		for(int i=1;i<aM.length;i++)
			for(int j=1;j<aM[1].length;j++){
				double sum=0;
				for(int k=1;k<aT[1].length;k++)
					sum+=aT[i][k]*a[k][j];
				aM[i][j]=sum;
			}
		
		for(int i=1;i<ABM.length;i++){
			double sum=0;
			for(int k=1;k<aT[1].length;k++)
				sum+=aT[i][k]*b[k];
			ABM[i]=sum;
		}
		
		LinearEquationGroup eqGroup=new LinearEquationGroup(n,aM,ABM);
		
		double[] solution=eqGroup.getSolution();
		
		for(int i=1;i<solution.length;i++)
			output+="x"+i+"="+fmt.format(solution[i])+"\n";
		
		return output;
	}

}