代码笔记 | java实现高斯赛德尔算法解线性方程组

package linear_equation;

import java.util.Scanner;

/*使用高斯赛德尔迭代法求解线性方程组*/
public class Gauss_Seidel_Iterate {
/*求下三角*/
private static float[][] find_lower(float data[][],int k){
int length=data.length;
float data2[][]=new float[length][length];
if(k>=0){
for(int i=0;i<=length-k-1;i++){
for(int j=0;j<=i+k;j++){
data2[i][j]=data[i][j];
}
}
for(int i=length-k;i<length;i++){
for(int j=0;j<length;j++){
data2[i][j]=data[i][j];
}
}
}
else{
for(int i=-k;i<length;i++){
for(int j=0;j<=i+k;j++){
data2[i][j]=data[i][j];
}
}
}
return data2;
}
/*求原矩阵的负*/
private static float[][] opposite_matrix(float[][] data){
int M=data.length;
int N=data[0].length;
float data_temp[][]=new float[M]
;
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
data_temp[i][j]=-data[i][j];
}
}
return data_temp;
}
/*原矩阵去掉第i+1行第j+1列后的剩余矩阵*/
private static float[][] get_complement(float[][] data, int i, int j) {

/* x和y为矩阵data的行数和列数 */
int x = data.length;
int y = data[0].length;

/* data2为所求剩余矩阵 */
float data2[][] = new float[x - 1][y - 1];
for (int k = 0; k < x - 1; k++) {
if (k < i) {
for (int kk = 0; kk < y - 1; kk++) {
if (kk < j) {
data2[k][kk] = data[k][kk];
} else {
data2[k][kk] = data[k][kk + 1];
}
}

} else {
for (int kk = 0; kk < y - 1; kk++) {
if (kk < j) {
data2[k][kk] = data[k + 1][kk];
} else {
data2[k][kk] = data[k + 1][kk + 1];
}
}
}
}
return data2;

}
/* 计算矩阵行列式 */
private static float cal_det(float[][] data) {
float ans=0;
/*若为2*2的矩阵可直接求值并返回*/
if(data[0].length==2){
ans=data[0][0]*data[1][1]-data[0][1]*data[1
4000
][0];
}
else{
for(int i=0;i<data[0].length;i++){
/*若矩阵不为2*2那么需求出矩阵第一行代数余子式的和*/
float[][] data_temp=get_complement(data, 0, i);
if(i%2==0){
/*递归*/
ans=ans+data[0][i]*cal_det(data_temp);
}
else{
ans=ans-data[0][i]*cal_det(data_temp);
}
}
}
return ans;

}

/*计算矩阵的伴随矩阵*/
private static float[][] ajoint(float[][] data) {
int M=data.length;
int N=data[0].length;
float data2[][]=new float[M]
;
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
if((i+j)%2==0){
data2[i][j]=cal_det(get_complement(data, i, j));
}
else{
data2[i][j]=-cal_det(get_complement(data, i, j));
}
}
}

return trans(data2);

}

/*转置矩阵*/
private static float [][]trans(float[][] data){
int i=data.length;
int j=data[0].length;
float[][] data2=new float[j][i];
for(int k2=0;k2<j;k2++){
for(int k1=0;k1<i;k1++){
data2[k2][k1]=data[k1][k2];
}
}

/*将矩阵转置便可得到伴随矩阵*/
return data2;

}

/*求矩阵的逆，输入参数为原矩阵*/
private static float[][] inv(float [][] data){
int M=data.length;
int N=data[0].length;
float data2[][]=new float[M]
;
float det_val=cal_det(data);
data2=ajoint(data);
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
data2[i][j]=data2[i][j]/det_val;
}
}

return data2;
}
/*矩阵加法*/
private static float[][] matrix_add(float[][] data1,float[][] data2){
int M=data1.length;
int N=data1[0].length;
float data[][]=new float[M]
;
for(int i=0;i<M;i++){
for(int j=0;j<N;j++){
data[i][j]=data1[i][j]+data2[i][j];
}
}
return data;
}
/*矩阵相乘*/
private static float[][] multiply(float[][] data1,float[][] data2){
int M=data1.length;
int N=data1[0].length;
int K=data2[0].length;
float[][] data3=new float[M][K];
for(int i=0;i<M;i++){
for(int j=0;j<K;j++){
for(int k=0;k<N;k++){
data3[i][j]+=data1[i][k]*data2[k][j];
}
}
}
return data3;
}
/*输入参数为原矩阵和一个整数,该整数代表从对角线往上或往下平移的元素个数*/
private static float[][] find_upper(float[][] data,int k){
int length=data.length;
int M=length-k;
float[][] data2=new float[length][length];
if(k>=0){
for(int i=0;i<M;i++){
for(int j=k;j<length;j++){
data2[i][j]=data[i][j];
}
k+=1;
}
}
else {
for(int i=0;i<-k;i++){
for(int j=0;j<length;j++){
data2[i][j]=data[i][j];
}
}
for(int i=-k;i<length;i++){
for(int j=i+k;j<length;j++){
data2[i][j]=data[i][j];
}
}
}
return data2;
}
/*m*n矩阵与n维向量的乘法*/
private static float[] multiply2(float[][] data1,float[] data2){
int M=data1.length;
int N=data1[0].length;
float[] data3=new float[M];
for(int k=0;k<M;k++){
for(int j=0;j<N;j++){
data3[k]+=data1[k][j]*data2[j];
}
}
return data3;
}
/*向量加法*/
private static float[] matrix_add2(float[] data1,float[] data2){
int M=data1.length;
float data[]=new float[M];
for(int i=0;i<M;i++){
data[i]=data1[i]+data2[i];
}
return data;
}
/*求两向量之差的二范数（用于检验误差）*/
private static double cal_error(float[] X1,float[] X2){
int M=X1.length;
double temp=0;
for(int i=0;i<M;i++){
temp+=Math.pow((X1[i]-X2[i]),2);
}
temp=Math.sqrt(temp);
return temp;
}
/*求矩阵的对角矩阵*/
private static float[][] find_diagnal(float A[][]) {
int m = A.length;
int n = A[0].length;
float B[][] = new float[m]
;
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
if (i == j) {
B[i][j] = A[i][j];
}
}
}
return B;

}
/*高斯赛德尔迭代法*/
private static float[] Gauss_Seidel_method(float[][] A,float[] B,float[] X){
float D[][]=find_diagnal(A);
float L[][]=find_lower(A, -1);
float U[][]=find_upper(A, 1);
float temp1[][]=inv(matrix_add(D, L));
float temp2[][]=opposite_matrix(temp1);
float B0[][]=multiply(temp2, U);
float F[]=multiply2(temp1, B);

return matrix_add2(multiply2(B0, X), F);

}

public static void main(String[] args) {
System.out.println("输入系数矩阵的行和列数：");
Scanner scan=new Scanner(System.in);
int M=scan.nextInt();
System.out.println("输入方程组右侧方程值的维度：");
int K=scan.nextInt();
if(M!=K){
System.out.println("方程组个数和未知数个数不等！");
System.exit(0);
}

System.out.println("输入系数矩阵：");
float[][] A=new float[M][M];
for(int i=0;i<M;i++){
for(int j=0;j<M;j++){
A[i][j]=scan.nextFloat();
}
}

System.out.println("输入值向量");
float[] B=new float[M];
for(int i=0;i<M;i++){
B[i]=scan.nextFloat();
}

System.out.println("输入初始迭代向量：");
float[] X=new float[M];
for(int i=0;i<M;i++){
X[i]=scan.nextFloat();
}

System.out.println("输入误差限：");
float er=scan.nextFloat();
float temp[]=new float[M];
while(cal_error((temp=Gauss_Seidel_method(A, B, X)), X)>=er){
X=temp;
}
//      while(cal_error((temp=Gauss_Seidel_method(A, B, X)), X)>=er){
//          X=temp;
//
//      }
X=temp;
System.out.println("高斯赛德尔计算得到的解向量为：");
for(int i=0;i<M;i++){
System.out.println(X[i]+" ");
}
System.out.println();

}
}