用C语言实现线性方程组的数值解法

简介：
数值求解方法主要有以下两条途径：
直接法：利用Gauss消元或矩阵分解，通过有限次运算可求出精确解。
迭代法：构造迭代格式，产生迭代序列，通过无限次迭代过程求解。有限次截断得近似解。
其中，直接解法可采用Gauss消元法、Doolittle分解法（三角分解法、LU分解法）、Crout分解法（与Doolittle类似）、Cholesky分解法（平方根法）。

标注：
此笔记中的四种消元方法已用函数封装，矩阵A和向量b数据的读取可以通过读文件操作，从文件A.txt与文件b.txt中进行读取（需自行创建），也可以通过数据输入的方式进行矩阵赋值。同时通过引用time.h头文件实现对四大算法的执行速度的测试。

程序代码：

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

//非线性方程的直接解法
void gauss(float*A,float*x,float*b,int n);
void Lu_doolittle(float*A,float*x,float*b,int n);
void Lu_crout(float*A,float*x,float*b,int n);
void cholesky(float*A,float*x,float*b,int n);

//输出数组元素
void printArray(float *a,int m,int n)
{
    for(int i=0; i<m; i++)
    {
        for(int j=0; j<n; j++)
        {
            printf("  %.2f",*(a+i*n+j));
        }
        printf("\n");
    }
}

//初始化数组，赋值为0
void initArray(float *a,int m,int n)
{
    for(int i=0; i<m; i++)
    {
        for(int j=0; j<n; j++)
        {
            *(a+i*n+j)=0;
        }
    }
}

void menu()
{
    printf("菜单：\n");
    printf("1. 读取矩阵A与矩阵b\n");
    printf("2. 输入矩阵A与矩阵b\n");
    printf("3. 使用GAUSS法解线性方程组(注：使用此方法后A内的值将改变)\n");
    printf("4. 使用Doolittle法解线性方程组\n");
    printf("5. 使用Crout法解线性方程组\n");
    printf("6. 使用Cholesky法解线性方程组\n");
    printf("0. 退出\n");
}

void InputMessage(float *a,int m,int n)
{
    for(int i=0; i<m; i++)
    {
        for(int j=0; j<n; j++)
        {
            scanf("%f",&a[i*n+j]);
        }
    }
}

int main()
{
    int n=0,i;
    float *A;//存放系数矩阵A
    float *b;//存放向量b
    float *x;//存放解向量x
    clock_t start, finish;//记录求解方程的开始时间和结束时间
    double duration;//保存求解方程所用的总时间

    //读取文件A与b
    FILE *fp1=fopen("A.txt","r");
    FILE *fp2=fopen("b.txt","r");

    if(fp1==NULL||fp2==NULL)
        return 0;

    menu();
    int choose;
    while(1)
    {
        printf("请输入你的选择：");
        scanf("%d",&choose);
        start=clock();//计时开始
        switch(choose)
        {
        case 1:
            printf("n元线性代数方程的n的值为：");
            scanf("%d",&n);
            A=(float*)malloc(n*n*sizeof(float));
            b=(float*)malloc(n*sizeof(float));
            x=(float*)malloc(n*sizeof(float));
            for(i=0; i<n*n; i++)
            {
                fscanf(fp1,"%f",A+i);//以浮点数的形式将文件中的数据取出赋值给矩阵A
            }
            for(i=0; i<n; i++)
            {
                fscanf(fp2,"%f",b+i);//以浮点数的形式将文件中的数据取出赋值给向量b
            }
            printf("输入读入成功！\n");
            printf("A=\n");
            printArray(A,n,n);
            printf("b=\n");
            printArray(b,n,1);
            break;
        case 2:
            printf("n元线性代数方程的n的值为：");
            scanf("%d",&n);
            A=(float*)malloc(n*n*sizeof(float));
            b=(float*)malloc(n*sizeof(float));
            x=(float*)malloc(n*sizeof(float));
            printf("输入矩阵A：");
            InputMessage(A,n,n);
            printf("输入向量b:");
            InputMessage(b,n,1);
            break;
        case 3:
            gauss(A,x,b,n);
            break;
        case 4:
            Lu_doolittle(A,x,b,n);
            break;
        case 5:
            Lu_crout(A,x,b,n);
            break;
        case 6:
            cholesky(A,x,b,n);
            break;
        case 0:
            exit(0);
            break;
        }
        if((choose!=2)&&(choose!=1))
        {
            finish=clock();//计时结束
            duration=(double)(finish - start);
            printf("用时%.2f ms\n",duration);
        }
    }
    //释放资源和关闭文件
    free(A);
    free(b);
    free(x);
    fclose(fp1);
    fclose(fp2);
    return 0;
}

//高斯消元法
void gauss(float*A,float*x,float*b,int n)
{
    float *t=(float*)malloc(n*n*sizeof(float));
    //矩阵化简
    for(int k=0; k<n; k++)
    {
        for(int i=k+1; i<n; i++)
        {
            t[i*n+k]=A[i*n+k]/A[k*n+k];
            b[i]=b[i]-t[i*n+k]*b[k];
            for(int j=0; j<n; j++)
            {
                A[i*n+j]=A[i*n+j]-t[i*n+k]*A[k*n+j];
            }
        }
    }
    //求解解向量x
    for(int i=n-1; i>=0; i--)
    {
        x[i]=b[i];
        for(int j=i+1; j<n; j++)
            x[i]-=(A[i*n+j]*x[j]);
        x[i]/=A[i*n+i];
    }
    printf("A=\n");
    printArray(A,n,n);
    printf("x=\n");
    printArray(x,n,1);
}
//LU分解法 杜里特尔
void Lu_doolittle(float*A,float*x,float*b,int n)
{
    float *L=(float *)malloc(n*n*sizeof(float));
    float *U=(float *)malloc(n*n*sizeof(float));
    float *y=(float *)malloc(n*sizeof(float));
    initArray(L,n,n);
    initArray(U,n,n);
    //初始化L的对角线上的值
    for(int i=0; i<n; i++)
        L[i*n+i]=1;
    for(int k=0; k<n; k++)
    {
        //计算U的第一行到第k行元素
        for(int j=k; j<n; j++)
        {
            U[k*n+j]=A[k*n+j];
            for(int i=0; i<=k-1; i++)
                U[k*n+j]-=(L[k*n+i]*U[i*n+j]);
        }
        //计算L中第一列到底第k列的元素
        for(int i=k+1; i<n; i++)
        {
            L[i*n+k]=A[i*n+k];
            for(int j=0; j<=k-1; j++)
                L[i*n+k]-=(L[i*n+j]*U[j*n+k]);
            L[i*n+k]/=U[k*n+k];
        }
    }
    //解方程组Ly=b
    for(int i=0; i<n; i++)
    {
        y[i]=b[i];
        for(int j=0; j<=i-1; j++)
            y[i]-=(L[i*n+j]*y[j]);
    }
    //解方程组Ux=y
    for(int i=n-1; i>=0; i--)
    {
        x[i]=y[i];
        for(int j=i+1; j<n; j++)
            x[i]-=(U[i*n+j]*x[j]);
        x[i]/=U[i*n+i];
    }
    printf("L=\n");
    printArray(L,n,n);
    printf("U=\n");
    printArray(U,n,n);
    printf("x=\n");
    printArray(x,n,1);
    free(L);
    free(U);
    free(y);
}
//LU分解法 克劳勃
void Lu_crout(float*A,float*x,float*b,int n)
{
    float *L=(float *)malloc(n*n*sizeof(float));
    float *U=(float *)malloc(n*n*sizeof(float));
    float *y=(float *)malloc(n*sizeof(float));
    initArray(L,n,n);
    initArray(U,n,n);
    //初始化U的对角线上的值
    for(int i=0; i<n; i++)
        U[i*n+i]=1;
    for(int k=0; k<n; k++)
    {
        //计算L的第一列到底第k列的元素
        for(int i=k; i<n; i++)
        {
            L[i*n+k]=A[i*n+k];
            for(int m=0; m<=k-1; m++)
                L[i*n+k]-=(L[i*n+m]*U[m*n+k]);
        }
        //计算U中第一行到第k行元素
        for(int j=k+1; j<n; j++)
        {
            U[k*n+j]=A[k*n+j];
            for(int m=0; m<=k-1; m++)
                U[k*n+j]-=(L[k*n+m]*U[m*n+j]);
            U[k*n+j]/=L[k*n+k];
        }
    }
    //解方程组Ly=b
    for(int i=0; i<n; i++)
    {
        y[i]=b[i];
        for(int j=0; j<=i-1; j++)
            y[i]-=(L[i*n+j]*y[j]);
        y[i]/=L[i*n+i];
    }
    //解方程组Ux=y
    for(int i=n-1; i>=0; i--)
    {
        x[i]=y[i];
        for(int j=i+1; j<n; j++)
            x[i]-=(U[i*n+j]*x[j]);
    }
    printf("L=\n");
    printArray(L,n,n);
    printf("U=\n");
    printArray(U,n,n);
    printf("x=\n");
    printArray(x,n,1);
    free(L);
    free(U);
    free(y);
}

//平方差法（要求矩阵A为对称正定矩阵）
void cholesky(float*A,float*x,float*b,int n)
{
    float *L=(float *)malloc(n*n*sizeof(float));
    float *U=(float *)malloc(n*n*sizeof(float));
    float *y=(float *)malloc(n*sizeof(float));
    initArray(L,n,n);
    initArray(U,n,n);
    L[0]=sqrt(A[0]);
    //求出矩阵L
    for(int k=0; k<n; k++)
    {
        for(int j=k; j<n; j++)
        {
            L[j*n+k]=A[k*n+j];
            for(int m=0; m<=k-1; m++)
            {
                L[j*n+k]-=L[k*n+m]*L[j*n+m];
            }
            if(k==j)
            {
                L[j*n+k]=sqrt(L[j*n+k]);
            }
            else
            {
                L[j*n+k]/=L[k*n+k];
            }
        }
    }
    //求出L的转置U
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
        {
            U[i*n+j]=L[j*n+i];
        }
    }
    //解方程组Ly=b
    for(int i=0; i<n; i++)
    {
        y[i]=b[i];
        for(int j=0; j<=i-1; j++)
            y[i]-=(L[i*n+j]*y[j]);
        y[i]/=L[i*n+i];
    }
    //解方程组Ux=y
    for(int i=n-1; i>=0; i--)
    {
        x[i]=y[i];
        for(int j=i+1; j<n; j++)
            x[i]-=(U[i*n+j]*x[j]);
        x[i]/=U[i*n+i];
    }
    printf("L=\n");
    printArray(L,n,n);
    printf("U=\n");
    printArray(U,n,n);
    printf("x=\n");
    printArray(x,n,1);
    free(L);
    free(U);
    free(y);

}

运行截图：