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][N];

         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 ][ 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][N];

         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][N];

         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][N];

         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][n];

         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();

 

     }

 

}

原文：http://www.oschina.net/code/snippet_574827_39084