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