C语言实现列主元消去法

372 阅读2分钟

​本文已参与「新人创作礼」活动,一起开启掘金创作之路。

 首先说明一下列主元消去法的内容,以下是360搜索到的,应该可以看懂

方法说明(以4阶为例):

第1步消元--在增广矩阵(A,b)第一列中找到绝对值最大的元素,将其所在

行与第一行交换,再对(A,b)做初等行变换使原方程组转化为如下形式:

第2步消元--在增广矩阵(A,b)中的第二列中(从第二行开始)找到绝对值

最大的元素,将其所在行与第二行交换,再对(A,b)做初等行变换使原方程组转化为:

第3步消元--在增广矩阵(A,b)中的第三列中(从第三行开始)找到绝对值

最大的元素,将其所在行与第三行交换,再对(A,b)做初等行变换使原方程组转化为:

按x4 → x3→ x2→ x1 的顺序回代求解出方程组的解。

然后C语言实现   参考了一位大神的博客 :   列选主元Guass消元法求解方程组+c语言_再__努力1点的博客-CSDN博客

//获取主元所在行
int getMaxinumLaber( double a[N][M], int k)
{
    int i;
    int laber=k;
    double maxinum=0;

    //列固定,行变,找最大值
    for(i=k; i<N; i++)
        if(maxinum<fabs(a[i][k]))
        {
            maxinum = fabs(a[i][k]);
            laber = i;
        }
    return laber;
}

//列主元高斯消去法
void gaussian_elimination( double a[N][M],double *x)
{
   int i, j, k;
    int laber;
    double temp;
   // double **a = a_mai;
    //double *x=x_mai;

    //高斯消元法
    for(k=0; k<N; k++)
    {
        laber = getMaxinumLaber(a, k);


        //主元交换,就是laber行和第k行交换
        if(laber!= k)
        {
          for(i=0; i<M; i++)
          {
              temp = a[k][i];
              a[k][i] = a[laber][i];
              a[laber][i] = temp;
          }
        }


        //消元
        for(i=k+1; i<N; i++)
        {

            if(a[k][k]==0)

            break;


            temp = a[i][k]/a[k][k];


            for(j=k; j<M; j++)
            a[i][j] = a[k][j]*temp-a[i][j];
        }
    }


    //回代求解x
    double sum;
    for(i=N-1; i>=0; i--)
    {
        sum = 0;
        for(j=i+1; j<N; j++)
        {
          sum += a[i][j]*x[j];
        }
        x[i] = (a[i][M-1]-sum)/a[i][i];
    }

    printf("\n");
    printf("变换增广系数矩阵\n");
    for(i=0; i<N; i++)
    {
        for(j=0; j<M; j++)
        printf("%f   ", a[i][j]);
        printf("\n");
    }


    printf("\n");
    for(i=0; i<N; i++)
    printf("%f   ", x[i]);
    printf("\n");

}

然后讲一下如何用:

a[N][M]其实就是我们线性代数中的增广系数矩阵,x[N]就是待求解的未知数向量;

应用是结合最小二乘法进行曲线拟合,首先通过最小二乘法求得增广系数矩阵,然后就是求拟合函数y =ax+b中的系数,b就是x[0],a就是x[1]   .

如果不了解最小二乘法的话就直接找个方程试试。

int main(void)
{
    double x_know[5] ={1,2,3,4,5},y_know[5] ={7,11,17,27,40},x[N]={0};
    int n =5,i=0,j=0;
    double A[N][M]={0};

    for(i=0;i<n;i++)
    {

      A[0][0]+=1;
      A[0][1]+=x_know[i];
      A[0][2]+=y_know[i];
      A[1][0]=A[0][1];
      A[1][1]+=x_know[i]*x_know[i];
      A[1][2]+=x_know[i]*y_know[i];
    }

    printf("\n");
    printf("原增广系数矩阵\n");
    for(i=0; i<N; i++)
    {
        for(j=0; j<M; j++)
        printf("%f   ", A[i][j]);
        printf("\n");
    }
    gaussian_elimination(A,x);

    return 0;
}