C语言 高斯-若尔当消元法

高斯若尔当方便解N元方程：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

float a[3][4]={{2,1,-1,8},{-3,-1,2,-11},{-2,1,2,-3}};
int rows=3,cols=4;


void print_matrix()
{//打印矩阵
    int i,j;
    int m=rows,n=cols;
    for(i=0;i<m;i++)
    {
        for(j=0;j<n;j++){
            printf("%10.3f\t",a[i][j]);
        }
        printf("\n");
    }
}
void swap_row(row1,row2)
{//交换行列
    int j=0,n=cols;
    float temp;
    for(;j<n;j++){
        temp=a[row1][j];
        a[row1][j]=a[row2][j];
        a[row2][j]=temp;
    }
}
void xiao_zhuyuan(int row,int curj,float val)
{//消主元，使之为1
    printf("行%d 除以 %f\n",row,val);
    int j=0,n=cols;
    for(;j<n;j++)
    {
        a[row][j]=a[row][j]/val;
    }
}

void gauss(){
    int i=0,j=0; //行和列
    int m=rows,n=cols;
    while(i<m&&j<n)
    {
            int k,maxi;
            float zhuyuan;
            k=i+1;

            if(fabs(a[k][j])>fabs(a[i][j]) ){

                swap_row(k,i); //按绝对值大小排序
                printf("行%d与行%d交换位置\n",k,i);
                print_matrix();
            }


            //消除非1的主元;
            zhuyuan=a[i][j];
            if(zhuyuan!=1){

               xiao_zhuyuan(i,j,zhuyuan);
                print_matrix();
            }//if



           if(a[k][j]!=0){
                int temp=0;
                int u,v;

                for(v=i;v<m-1;v++){
                float per=a[v+1][j]/a[i][j];//系数

                    printf("行%d减 (%f倍)行%d\n",v+1,per,i);
                    for(u=0;u<n;u++){
                    a[v+1][u]=a[v+1][u]-per*a[i][u];

                     }
                      print_matrix();

                }
                 j++;
            }
             i++;


    }//while
}//gauss

int main()
{
    //高斯消元法的算法复杂度是O(n3);这就是说，如果系数矩阵的是n × n，那么高斯消元法所需要的计算量大约与n3成比例

    /*  2x+y-z=8
     *  -3x-y+2z=-11
     *  -2x+y+2z=-3
     * 因此，增广举证为：
     *  2   1   -1      8
     * -3   -1  2       -11
     * -2   1   2       -3
     */

     int i,j;
     printf("原始增广矩阵\n");
     print_matrix();
     gauss();

     //printf("消元后的矩阵\n");
    // print_matrix();

}

不过目前只做到 矩阵梯队的形式。

继续。。。。。

——————————————————————————————————————————————————————————————————

上个例子中存在一些问题：

1.在交换行的时候，没有全行比较，只比较了相邻2行。

2.在显示第几行的时候，用的是数组的下标，故少了1.

3.没有对浮点数的显示优化，现在用的是"%.3g"，这样能减少小数部分的长度.

4.没有用到 jordan的方法直接求解。

下面是完整例子，不过也不是最佳的，因为里面的变量看起来有点乱。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

float a[3][4]={{2,1,-1,8},{-3,-1,2,-11},{-2,1,2,-3}};
//float a[4][5]={{2,1,-1,8,1},{-3,-1,2,-11,1},{-2,1,2,-3,1},{5,7,8,0,1}};
int rows=3,cols=4;


void print_matrix()
{//打印矩阵
    int i,j;
    int m=rows,n=cols;
    for(i=0;i<m;i++)
    {
        for(j=0;j<n;j++){
            printf("%10.3g",a[i][j]);
        }
        printf("\n");
    }
     printf("\n");
}
void swap_row(row1,row2)
{//交换行列
    int j=0,n=cols;
    float temp;
    for(;j<n;j++){
        temp=a[row1][j];
        a[row1][j]=a[row2][j];
        a[row2][j]=temp;
    }
}
void xiao_zhuyuan(int row,int curj,float val)
{//消主元，使之为1
    printf("行%2d 消除 %.3g\n",row+1,val);
    int j=0,n=cols;
    for(;j<n;j++)
    {
        a[row][j]=a[row][j]/val;
    }
}

void gauss(){
    printf("\n高斯消元开始:\n");
    int i=0,j=0; //行和列
    int m=rows,n=cols;
    while(i<m&&j<n)
    {
            int k,d,maxi=i,maxval;
            float zhuyuan;
            k=i+1;
            maxval=fabs(a[i][j]);
            for(d=k;d<n;d++){
                 if(fabs(a[d][j])>maxval ){
                     maxval=fabs(a[d][j]);
                     maxi=d;
                }
            }
            if(maxi!=i){
               swap_row(maxi,i); //按绝对值大小排序
                printf("行%2d与行%2d交换位置\n",maxi+1,i+1);
                print_matrix();
            }


            //消除非1的主元;
            zhuyuan=a[i][j];
            if(zhuyuan!=1){
               xiao_zhuyuan(i,j,zhuyuan);
                print_matrix();
            }//if



           if(a[k][j]!=0){
                int temp=0;
                int u,v;

                for(v=i;v<m-1;v++){
                float per=a[v+1][j]/a[i][j];//系数

                    printf("行%2d = 行%2d 减去(%.3g )倍的行%2d\n",v+2,v+2,per,i+1);
                    //printf("Subtract (%.3f × row %d) from row %d\n",per,i+1,v+2);
                    for(u=0;u<n;u++){
                    a[v+1][u]=a[v+1][u]-per*a[i][u];
                     }
                    print_matrix();
                }//for
                 j++;
            }//if
             i++;
    }//while
}//gauss

void jordan()
{
printf("\n若尔当消元开始:\n");
 int i=0,j=0,k; //行和列
    int m=rows,n=cols;
    float xs;
    for(i=m-1;i>=0;i--){
        for(k=0;k<i;k++){
            xs=a[k][i];
            if(i!=k){
            printf("行%2d = 行%2d 减去(%.3g )倍的行%2d\n",k+1,k+1,xs,i+1);
            //printf("Subtract (%.3f × row %d) from row %d\n",xs,i+1,k+1);
              for(j=0;j<n;j++){
                a[k][j]=a[k][j]-xs*a[i][j];
                }
               print_matrix();
             }//if
       }
    }//for i

}//jordan
int main()
{

     int i,j;
     printf("原始增广矩阵\n");
     print_matrix();
     gauss(); //使用高斯消元
     jordan(); //使用若尔当消元
}

 

在测试这个程序的时候发现：如果是方阵话，若尔当消元后，变成类似一下的矩阵：

1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 1

附：

高斯消元法百科：http://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%B6%88%E5%8E%BB%E6%B3%95
高斯若尔当消元法百科：http://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF-%E8%8B%A5%E7%88%BE%E7%95%B6%E6%B6%88%E5%85%83%E6%B3%95
在线高斯-若尔当消元计算（英文）：http://www.idomaths.com/gauss_jordan.php